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ABSTRACT 

The extraction and upgrading of bitumen at _ surface 
mining projects require the construction of heat generating 
structures. In some cases these structures must be situated 
above oil sand strata. Oil sand expands upon heating with 
the amount of expansion dependent on the amount of drainage 
of water and bitumen from the heated zone. The expansion can 
cause large surface movements which can have adverse effects 
on the structures. 

Laboratory work consisted of the measurement of 
geotechnical properties of oil sand subjected to high 
temperatures. Detailed information was collected with regard 
to the influence of temperature on properties such as 
permeability, compressibility and volume change. Laboratory 
work also involved the preparation of high quality 
undisturbed oil sand test samples from frozen core. 

Analytical studies involved Suta lazing! \antormation 
derived from the laboratory testing program to investigate 
theoretical methods for predicting volume change and pore 
fluid development of oil sand subjected to increases in 
temperature. Numerical methods were employed to analyze the 
rate at which heat flows through an oil sand mass and the 
rate at which excess pore fluid pressure resulting from an 
increase in temperature dissipates. A method was developed 
to couple the above processes (heat transfer and fluid 


flow). 
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Undrained heating of 011 sand to 200 °C will result in 
a volume increase of about 6 percent. Most of this increase, 
over 5 percent, iS due to the pore fluids, water and 
bitumen, with the bitumen expanding slightly more than 
water. The remainder of the volume increase is due to _ the 
thermal expansion of the sand grains. The aenser 
Pigehechanine Structure of in situ oil sand will show a 
negligible amount of thermal volume change by itself. 

The analytical model presented in this study 
investigated the effect of heat flow from a 80 m diameter 
Storage tank on underlying oil sand strata. The tank was 
assumed to rest on a Sue heaceal crushed rock and sand 
foundation pad, 2 m in thickness, and held. at a constant 
semperature ‘of 175°C; 

The results of the theoretical analysis showed that 
although the rate of pore fluid flow was slower than that of 
heat flow, it was sufficient to provide a significant amount 
of drainage. The pore fluid pressure developed during 
heating of the oil sand strata was found to drain fast 
enough to negate large volume changes associated with the 
undrained condition. A maximum differential heave of 14 cm 
was observed after 30 years of heating. The amount of 
differential heave predicted from the analysis appears 
manageable when incorporated into a structural design of a 


- storage tank. 
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1. INTRODUCTION 


1.1 Statement of Problem 

The purpose of this thesis is to examine the influence 
that temperature increases from heated foundations have on 
the foundation performance of oil sand and to determine 
whether unacceptable foundation movement could be caused to 
overlying structures. The aim of the study is also to 
determine what temperatures are allowable in the oil sand or 
what rate of temperature increase in the oil sand from 
heated foundations is allowable before alteration of the 
physical properties of the oil sand occurs. 

To analyze this problem, it 1S necessary to assess the 
rate at which heat generated by the structure propagates 
through the underlying oil sand strata and evaluate the 
amount of drainage which occurs. If drainage occurs at a 
rate greater than or equal to that of heat flow, the oil 
Sand mass will be fully drained and thermal volume changes 
are due to the sand grains and soil structure. If drainage 
occurs at a rate less than that of heat flow the oil sand 
mass will be partially drained and thermal volume changes 
are due to the sand grains, soil structure and pore fluids 
(bitumen and water). 

Lo tpdricorm ithe vtheat 6 flow Vtandlysis, * the -heat Bilow 
parameters are required. These include thermal conductivity, 
density and specific heat. Oil sand subjected to an increase 


in temperature will experience: an increase in volume; a 
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change in viscosity of the pore fluids; and an increase in 
pore fluid pressure. 

To determine the state of drainage it is necessary to 
assess: the distribution of pore pressure increase; 
permeability (which is dependent on viscosity); increase in 
volume of pore fluid; compressibility of pore fluid; and the 
change in volume of soil structure. 

The amount of heave and differential heave occurring 
beneath the heated structure are calculated from the volume 


changes taking place in the oil sand mass. 


142 Scope of Thesis 

Undisturbed onl sand Satisfies the fundamental 
requirements necessary for a good foundation material. Owing 
to its extremely dense interpenetrative structure, oil sand 
in situ has a high resistance to shearing failure 
(Chapter 2) and is virtually incompressible enabling it to 
Support masSive Structures with minimal settlements. To 
illustrate the latter, consider a storage tank 80 m in 
diameter containing 100,000 m*? of bitumen (density = 1.05 
Mg/m?) resting on a layer of oil sand 80 m in thickness. 
From the theory of elasticity (Lambe and Whitman, 1969) the 
stress at mid-layer (tank center) induced by the uniform 
nermalsstress FOre@b205 .kPaseti 05% MaYm° Stk 10020008%m" *% 
OXBA1/SGAe Nagel 4onagiled: <is FAAS kPaANC0F64*x5185HF “ASSiming@a 
coefficient of volume change of 5 x 10°’ kPa™' (Dusseault, 


1980) the settlement under the center of the tank would be 
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Oil sand subjected to temperature increases will expand 
with the interstitual fluids (water, bitumen and gases) 
contributing to the majority of the volume increase. As heat 
from the foundation propagates through the o11 sand, the 
thermal expansion will produce heaving beneath the 
foundation. The severity of the heaving is dependent on the 
EompoSsit@ionalmnature onthe cil escand Jocal. stratigraphy, 
foundation temperature and whether or not’ sufficient 
drainage of pore fluids can take place. 

The structure considered in the analyses is an 80m 
diameter bitumen storage tank maintained at a constant 
temperature of 175 °C (Figure 1.1). 

When examining the foundation performance. of the oil 
sand beneath a heated structure, the possibility of shear 
failure, lateral movement and settlement must also be 
considered. The possibility of shear failure is governed by 
the lateral confining pressure on the 011 sand beneath the 
Structure, the relative density of the oil sand and the size 
of the foundation. In the surface mineable area the minor 
principal stress is believed to be vertical from the surface 
to a depth of approximately 330 m and is one-third that of 
the horizontal stress to a depth of about 100 m (Dusseault, 
1977b)... The Mrelamivemidensitvyve.sot sine situ oll sand. is 
extremely high (chapter 3). Since the size of the foundation 
considered in this study is extremely large and the lateral 


confining pressure and relative density high, the 
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possibility of shear failure is extremely remote. These 
factors also negate the possibility of the oil sand and 
structure moving laterally. Results of this study will also 
show that settlement of the structure does not take place 
even under successive heating and cooling cycles. Heave and 
differential heave is the sole factor affecting the 
foundation performance of the oil sand. 

The magnitude of expansion that oil sand and its 
individual constituent components undergo with temperature 
increase was determined experimentally in the laboratory. 
Thermal volume change experiments were performed on 
undisturbed oil sand under completely undrained and drained 
conditions oe temperatures. ./of, 200 4.°C. -Bach component 
contributing to volume change; water, bitumen, sand grains, 
and soil structure were examined separately. The 
experimental values were compared to theoretical models for 
completely undrained and drained states. Experimental 
thermal constants derived from the testing of the oil sand 
constituent components were also used to examine the 
validity of existing formulations, theoretically derived for 
the prediction of pore pressure generation in heated soils. 

The analytical model developed in this study couples 
heat flow theory to fluid flow to assess the state of 
drainage in the oil sand strata below the heated foundation 
and calculate the change in pore pressure and the volume 
change (ormpheave)iq The eprepegation pofgrheatrefromenthe 


overlying structure through the oil sand Strata was 
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determined analytically utilizing numerical heat transfer 
techniques. With consideration to changes in permeability, 
coefficients of thermal volume change and compressibility, 
dynamic viscosity, density, thermal generation rate and 
specific storage coefficient of oil sand with temperature 
and preGsurer the amount of pore pressure drainage was also 
examined by numerical methods. 

The analytical model used in this study was developed 
for the two-dimensional case. Although the model is easily 
adapted to the three-dimensional case, results of this study 
will show that the increase in accuracy would not justify 


the large computational effort and expense. 


1.3 Organization of Thesis 

Published information concerning the characteristics of 
the Athabasca Oil Sands (including compositional nature and 
engineering properties), surface mining projects, and 
Studies relevant to heated foundations on oil sand is 
presented in Chapter 2. 

Chapter 3 contains a discussion of the analytical 
concepts utilized in this study. These concepts include: the 
influence of heat on oil sands (expansion and drainage); 
analysaseeotrihneat® Plowss analysis! of fluid flow; and* the 
coupling (Ot@heateand fluid flow: 

A description of the laboratory apparatus and test 
procedures employed in the experimental program are given in 


Chapter 4. The experimental results and discussion of 
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results are contained in Chapter 5. 

Numerical analysis simulating heat flow into an oil 
sand foundation from a hot bitumen storage tank is given in 
Chapter 6. Also in this chapter analytical techniques are 
used to predict the amount of drainage which may occur. The 
resulting volume change of the oil sand foundation from the 
numerical analyses is also given. 

Chapter 7 contains observations made during the study 
and conclusions drawn from the research. 

The engineering properties of the oil sand constituent 
components which were required as input for the numerical 
analyses are given in Appendix A. This appendix outlines the 
source of the data used to derive each property 
(experimental testing or published information) and how each 
property varies with temperature and pressure. 

Appendix B provides a detailed description of the 
apparatus used in the experimental testing program. Appendix 
C outlines the laboratory calibrations used and Appendix D 
provides the data reduction calculations. 

Appendix E contains listings of the available computer 


programs used in the study. 
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2. LITERATURE REVIEW 


2.1 The Athabasca Oil Sands 

The Athabasca Oil Sand Deposit is located in 
northeastern Alberta (Figure 2.1) underlying an area of 
approximately 32,000 square kilometers with estimated 
inplace reserves of 146.5 billion cubic meters (Outrim and 
Evans, 1978; Mossop, 1980). Athabasca, the largest of the 
four major oil sands See in Alberta, 1S unique in that 
a portion of the reserves (about 10 percent) have less’ than 
50 m of overburden (figure 1.1)and is exploitable by surface 
mining methods (Mossop, 1978). 

In contrast to the other Alberta oil sand deposits, the 
Athabasca reserves are limited for the most ‘part to one 
reservoir unit, the McMurray Formation. Mossop (1980) 
describes the McMurray Formation as a Lower Cretaceous unit 
consisting of uncemented quartz sand with interbedded 
Shales. In general, the entire McMurray Formation bears oil 
with the well sorted, shale free oil sand layers having 
porosities of up to 35 percent and oil saturations of 18 


weight percent (Carrigy and Kramers, 1974; Mossop, 1978). 


2.2 Characteristics of the Athabasca Oil Sands 

The vast majority of reserves of the Athabasca Deposit 
are found in the sands of the Lower Mannville McMurray 
Formation. The McMurray Formation in the area of concern 


averages between 40 and 60 m in thickness (Mossop, 1980). In 
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Figure 2.1 Oil Sands Deposits of Alberta 
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the area of the present mining operations, the overburden 
varies in thickness from 0 to about 35 m. 

Oil sand typical of the Athabasca deposit consists of 
approximately 95 percent quartz, 2 to 3 percent feldspar 
grains, 2 to 3 percent mica and clay minerals, and traces of 
other minerals (Mossop, 1978). The grains are primarily 
Subangular with moderate sphericity (Mossop, 1980). The high 
grade oil bearing sands in Athabasca are predominantly fine 
to medium grained and uniformly graded quartz sand 
(Dusseault, 1977). 

The sand grains are surrounded by a thin film of water, 
with the remaining pore space occupied by the oil 
(Figure 2.2). The water and the oil form continuous’ phases 
(Mossop, 1980). The majority of clay particles in the oil 
Sand matrix are contained in the water envelope usually 
attached directly to the sand grains. Substantial portions 
of gases (mostly methane and carbon dioxide) are dissolved 
within the liquid phase of the oil sands (Dusseault, 1977). 

Dusseault (1977; 1980) reports saturated bulk densities 
for uniformly graded, rich oil sand (the ore bearing zones 
of the Middle and Lower McMurray Formation) of 2.05 to 2.18 
Mg/m? and porosities of 28 to 36 percent. Mossop (1978) 
reports oil and water saturations of rich oil sands in the 
Athabasca Deposit of 18 and 2 weight percent, respectively. 
The grain size distribution for Athabasca oil sand is given 
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Figure 2.2 Composition of In Situ Oil Sand 
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Undisturbed oil sand exhibits extremely high shear 
strengths and dilatancy compared to normal dense sand of 
Similar minerology. Dusseault (1977) provides evidence from 
thin section and scanning electron microscope investigations 
which show oil sand, unlike typical dense sands, has a large 
number of concavo-convex and straight (or long) contacts. 
This interpenetrative structure developed during the time of 
burial in which diagentic processes occurred with 
dissolution and recrystall izataonieofs i quartzihatcegrain 
boundaries. 

Undisturbed oil sands exhibit curvilinear Mohr failure 
envelopes. The internal anghes ofifriction masa highsaat wlow 
normal stresses, reflecting high rates of dialation but as 
normal stresses increase dilatancy is reduced as more shear 
occurs through grains and surface asperities. The rate of 
eurvaturetin’ thes failure envelope decreases as the residual 
angle of shearing resistance is approached (33 to 36° 
(Dusseault, 1977). Dusseault (1977) expresses the failure 


envelopes for undisturbed oil sand as a power law functions. 


2.3 Surface Mining Projects 

Of the 11.8 billion cubic meters of reserves which lie 
within the surface mineable area of the Athabasca Deposit 
only 6.4 billion cubic meters are recoverable by current 
surface mining methods, resulting in 4.3 billion cubic 
meters of synthetic crude oil (Outrim and Evans, 1978). 


Figure 2.4 illustrates the surface mineable area of the 
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Athabasca oil sand deposit indicating the location of the 
two existing commercial operations (Suncor Inc. and Syncrude 
Canada Ltd.). A third surface mining project was proposed by 
the Alsands Consortium and was to be located 70 km north of 
Fort McMurray (Figure 2.4). 

Suncor began surface mining of oil sand with bucket 
wheel excavators in 1963 and is currently producing about 
BO000 cubic#* metersi#fot Fsynthetic ¥crude.o1)]) per day. The 
surface mining activity at the Syncrude site commenced in 
1977 with draglines and is currently producing approximately 
20,000 cubic meters of synthetic crude oil per day. The 
proposed Alsands project was planned to mine oil sand with 
draglines and was projected to produce approximately 22,000 


cubic meters of synthetic crude oil per day. 


2.4 Effect of Heat on O11 Sands 

Only recently, upon the construction of the Syncrude 
Canada Ltd. commercial production facility was the necessity 
Or erecting heated structures on’ orl “sands strata 
considered. As a result, little information exists in the 
published literature concerning the behavior of oil sands 
beneath heated foundations. 

Boga, et al. (1980) conducted a study to explore the 
StaDime ty Ofmcll sana exposed @uo excessive,, Nedting , beneath 
hot bitumen storage tanks. The investigation centered on the 
proposed installation of,.four additional hot storage tanks 


12 m above oil sand at the Syncrude plant site. The study 
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Figure 2.4 Surface Mineable Area of the Athabasca Oil Sands Deposit 
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waS intended to provide an optimum thermal foundation design 
for the hot bitumen storage tanks. 

Utilizing a finite element thermal model Boga, et al. 
(1980) investigated the propagation of heat from the tanks 
in the subsurface material. Triaxial testing of oil sand 
Samples under drained and undrained conditions to 
temperatures of 90 °C were also conducted. The o11 sand used 
in the testing program were obtained from block samples 
taken from the Syncrude mining operation. Since no 
percautions were taken to guard against expansion of the oil 
sand from in situ conditions, the samples could be assumed 
to be at a considerably lower density and thus be considered 
disturbed. The in situ horizontal stress in this study was 
assumed to be one-half the vertical in the Syncrude mine 
area. In the surface mineable area, however, the horizontal 
stress is believed to be three times the vertical stress to 
a depth of about 100 m (Dusseault, 1977b; section 1.2). Such 
conditions suggest one-dimensional testing (zero lateral 
Strain) is more representative of actual field conditions. 

The triaxial tests conducted in this study indicate 
that under fully drained conditions oil sand did not display 
anyinsoiteningrcomafakliures tometemperatures' of 90 °C. The 
triaxial tests conducted under fully undrained conditions, 
however, displayed a significant increase in pore fluid 
pressure at 60 °C and softened and subsequently failed at 70 
°C. It was concluded from these tests that under  undrained 


conditions the strength of oil sand is adversely affected at 
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temperatures exceeding approxémately 45IS8C. 

Boga, et al. (1980) noted in their investigation that 
considering the slow rate of heat flow into the oil sand 
foundation beneath the heated storage tank, some dissipation 
of pore fluid pressure was possible. However, since the 
amount of pore pressure dissipation from drainage could not 
be calculated, a maximum limiting temperature of 45 °C 
(based on undrained tests) for the oil sand beneath the tank 
was assumed for design purposes. This limiting temperature, 
combined with the thermal analysis, provided a foundation 
design for the proposed tanks. 

Harris and Sobkowicz (1977) developed a mathematical 
model to analyze the undrained behavior of oil sand to 
changes in temperature (and/or stress). The model has_ the 
capability of dealing with either one-dimensional 
compression (no lateral strain) or two-dimensional, plane 
Strain. The model does not, however, incorporate the effect 
of fluid pressure drainage which occurs during the time of 
heat propagation into the oil sand layer. As a result, the 
model may exaggerate the magnitude of potential heave which 
may occur as the temperature of the oil sand beneath the hot 
foundation increases with time. 

Byrne, et al. (1980) provide a procedure for estimating 
the undrained behavior of oil sand resulting from changes in 
temperatures and loads. The method involves use of the 
computer program OILSTRESS (developed by Byrne and Grigg 


[1980]). The program utilizes finite element techniques and 
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has the capability to deal with both plane strain and 
axisymmetric conditions. The program assesses volume changes 
occurring in the pore fluid and couples it with that of the 
sand matrix to maintain volumetric compatibility. An 
iterative procedure is used to obtain compatibility in the 
sand skeleton and fluid phases in terms of volumetric 
Strain. Gas laws are used to calculate changes in pore 
pressures. In the analysis the sand skeleton is considered 
non-linear elastic. The program also considers volume 
changes arising from shear dilatation: 

Although Byrne, et al. (1980) have not applied the 
computer program OILSTRESS directly to heated foundations on 
oil sand, it has been used successfully in predicting the 
benavioriof ‘one-dimensionalstmesstreliefvot oil sand” The 
procedure has also been utilized to predict the response of 
a deep shaft in oil sand. When applied to cylindrical 
openings in elastic and elastic-plastic materials the 
results obtained from the analysis compare favourably with 
closed form solutions. 

Charlwood, et al. (1980) modified the computer model 
put "cforthuibys'Byrne, et al.(4980) (which only considers 
temperature influence on the gas phase of oil sand) to 
include transient heat transfer and thermal expansion 
effects. The model is intended for use in mine assisted in 
Situ (MAIS) design applications. 

The computer model proposed by Charlwood et al. (1980) 


was not used for this study for several reasons. The model 
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assumes in itS pore pressure generation term for undrained 
conditions that water and bitumen are incompressible. It 
will be shown in this study that the compressibility of 
bitumen and water is similar to that of the soil Structure. 
The pore pressure generation term also does not include 
thermal expansion 5h pore fluids. 

The forces associated with temperature are calculated 
in this model based on a single value of thermal expansion 
coefficient and therefore fails to differentiate between 
drained and undrained thermal expansion. Therefore, if an 
undrained thermal expansion coefficient were used in the 
analysis, the stress distribution in the oil sand _ strata 
would only be correct if all elements remain fully undrained 
throughout the analysis. Should any of the elements reach a 
temperature in which the permeability was such to permit 
drainage of pore fluids, these elements would have to be 
manually altered to the correct drained volume in order to 
produce the proper stress distribution. In other words, the 
computer program does not account for mass transfer during 
Simultaneous thermal expansion and drainage because it 
cannot automatically differentiate between drained and 
undrained coefficients of thermal expansion. AS ae result, 
use of the program for the purposes required for this study 
would be extremely difficult. 

The geotechnical aspects of oil sand behavior due to 
temperature changes and the concept of heat-consolidation 


(the time rate of drainage during heating) are discussed by 
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Morgenstern (1981). 

A theoretical analysis for predicting volume changes in 
soils resulting from temperature variations and _ pore 
pressure changes due to undrained heating have been 
developed by Campanella and Mitchell (1968) and examined 
further by Mitchell (1976). These theoretical expressions 


are discussed in Chapter 3. 
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3. ANALYTICAL CONCEPTS 
3.1 Influence of Heat on Oil Sands 


3.1.1 Expansion 

A theoretical analysis for predicting volume changes in 
Saturated soils subjected to temperature changes was put 
forth by Campanella and Mitchell (1968). This analysis is 
extended herein to consider oil sand which contains two pore 
fluids (bitumen and water). The analysis presented here 
assumes the pressure of the pore fluids is sufficient to 
prevent gas exsolution and pore fluid vapourization or 
distillation. The influence of possible gas exsolution is 


discussed in the conclusions of Chapter 6. 


3.1.1.1 Expansion Under Fully Drained Conditions With no 
Effective Confining Pressure 

Under fully drained conditions the volume change of 
oil sand subjected to increases in temperature under no 
effective confining pressure arises from the expansion 
of the individual sand grains and structural changes in 
the sand skeleton. 

Campanella and Mitchell (1968) theorize that when a 
homogeneous sand with grains in particle to particle 
contact undergoes a temperature change, each mineral 
experiences the same volumetric strain and the entire 
Structure will Pee nence the same volumetric strain 


(a,4T). It is believed that subsequent volume changes 


Pa 


° 
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may arise from changes in intergranular forces due to 
temperature which entails relative movement or 
reorientation of particles to allow the sand mass to 
Support the same effective stress. Campanella and 
Mitchell (1968) denote this change in volume by (AV,,),.. 
As a result the net volume change of an oil sand mass 
subjected to a change in temperature under fully drained 


conditions is: 


(AVm) = a2 VAT (AV. +) oa 
where: 
a, = volumetric coefficient of thermal expansion of the 


sand grains 
Vm = volume of oil sand mass 


AT = change in temperature 


The structural volume change of the sand structure due 


to a change in temperature may be expressed as: 


(AVG ae aT ee 
where: 
@s;+ = coefficient of structural volume change of sand 


structure due to change in temperature 


4 7 hae : aS r a: a : vote —o ia ie > 


it 5 th) Ab \ ; : ; 
f ret } ; a 


ao sub 205702 ssluneapasaat i aavheits morta ‘ed 


Hie iv 


~~ tmp vce : svitetor | alieine doiiy. ore 


| od Be Bi bnse ‘aiid woke oo nokiateaed tS ves ceon 
ria . ati enaqme? De qashtz ev iiverte shins | ont | tae 1 
“ris ; wal ¥a amu Lov et spnieit's atid 92 rama (aye i Lae ei 
| oe bras fio: «te FO Sef ais pap low feet art aiuese 
beni arb gl ‘ ui seBau. seoeeS damnect no sored — od Reds she 
, | +g. 
| fet ano b oie 
5 ; = yy 
eid 29 fodermgee Luger ert ?O thei ott gen serteaed ea *” 
| 226m hbase [toa ia smtey = “ 
STU te) some oc? spneis, = 
e 7 | | | 
et $243 we Snee: ens Ro Spee SM ary) leuusau ste edt « 
¢ > 
: tes Besteides sd vem sou mass ni “sprite & +) 
a | 
Po fi : hd Gs 3 | al 
: S.£ : a pe ) TAY fue rahe 


a re 
See 7 ‘i 


23 


Combining equations 3.1 and 3.2 gives: 


CAVig) = Calg et a. «VV gAT 333 


3.1.1.2 Expansion Under Completely Undrained Conditions 
With Constant Effective Stress 

Under completely undrained conditions the volume 
change of oil sand subjected to increases in 
temperatures under a constant effective stress arises 
from the expansion of the individual sand grains and 
pore fluids as well as structural changes in the_- sand 
skeleton. 

The volume change of pore fluids experiencing a 


change in temperature can be stated mathematically as: 


(AVE) = CV CAT Rowe G 3.4 


wheres: 


a,= coefficient of thermal volume change of water 
Vw= volume of pore water 

AT = change in temperature 

a,= coefficient of thermal volume change of bitumen 


Vp= volume of pore bitumen 


The total volume change of an oil sand mass 
Subjected to a change in temperature under completely 


undrained conditions (effective stress held constant) 
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is: 

os d 
(AVm) = (AV) (AVm"), + BBs) 
wheres 
(AVm’?) = net volume change of an oil sand mass subjected 
to a change in temperature under fully drained 
conditions 


Combining equations 3.3, 3.4, and 3.5 gives: 
(AN) 20 aM, ATs to, Vea ety (a) + Clee) VA 306 


3.1.2 Volume of Pore Fluid Expelled under Fully Drained 
Conditions with No Effective Confining Pressure 
The volume of pore fluid Anes (under no confining 
stress) from a fully saturated oil sand sample experiencing 


a change in temperature can be expressed as: 


(AVe) = (AVE) + (AV) > ( AVm ler aha 


where: 


(AV, ),-= volume change of the sand grains due to change in 


temperature 
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The volume change of the sand grains due to change in 


temperature can be expressed mathematically as: 


(AViO =wial Viv .OT 3.8 


Combining equations 3.3, 3.4, 3.6 and 3.7 gives: 


(AV 4 ae on Ve oTe + a. Vy, ADS te ioe VL AT  ehVeATorcoa.ce Vi, AT Sao 
ors 

(AVa) = ae Va dtote af Votes — a, AT(Ve,-Vi,) Wak VAT 3.10 
where: 

Via -— V, = volume of the voids 


Put simply into words, the amount of fluid expelled 
from an oil sand sample during a drained thermal expansion 
test is the sum of the volume change of the fluids less the 
change in void volume and structural volume change. If the 
structure expands [AV.,] is positive and would be subtracted 
Or contractungly aut e'the: stnuetures contracts, [AV,.J,7as 


negative and would be added. 
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3.1.3 Heat Consolidation 

The amount of expansion an oil sand mass experiences 
due to temperature changes depends on the amount of pore 
fluid which is permitted to drain. To evaluate the quantity 
of fluid flowing from oil sand subjected to temperature 
increases it is necessary to assess the relative rates of 
pore pressure generation and dissipation. 

Campanella and Mitchell (1968) provide a theoretical 
analysis for predicting pore water pressure changes in 
Saturated soils subjected to temperature changes. This 
analysis is extended here to consider oil sand which 
contains two pore fluids (bitumen and water). Gas exsolution 
and pore ..flwid. vapourization is not considered in -the 


analysis. 


3.1.3.1 Increase in Pore Fluid Pressure 

Under undrained conditions the combined volume 
changes of the oil sand constituents (bitumen, water and 
Sand grains) resulting from changes in pressure and 
temperature must balance the combined changes of the 
entire oil sand mass resulting from changes in pressure 


and temperature, i.e.: 


(AV, ) 


rt (AV 


yt (AV. et (AVE), * (AV sgt (AVE y= (BV art 


(AVimn) 4, Beni 
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AV, = volume change of pore bitumen 
AVw= volume change of pore water 


AV, 


volume change of sand grains 


AVE volume change of oil sand mass 
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The volume changes of the oil sand constituents due 


to pressure changes can be written as: 
(AV), = m, V, Au 

(AV ap= MyVy Au 

(AV, ao= m,V,Au + m,'V,Ao! 

where: 

m= compressibility of bitumen 


Au = change in pore fluid pressure 


my= compressibility of water 


m, = compressibility of sand grains under an 
pressure 
m,’ = compressibility of sand grains 


subjected to concentrated loadings 


Ao’. = change,an intergranular (effective) 


round 


are 


The volume change of the total oil sand mass due to 


pressure changes can be written as: 
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(AV y= My Vmdo' 3.15 
where: 
m,= compressibility of oil sand mass 


The change in volume of pore water, pore bitumen 


and sand grains due to temperature change are expressed 


ass: 
(AV, ) = awVy AT ore 
(AV, ) = a, Vp AT Jihad 


Substituting equation 3.8 and equations She 
through 3.17 into the governing equation for undrained 


conditions (equation 3.11) gives: 


ae Mp YATULt ag yaTiona dv saTat mp Vp Au Sane AUG + meV oe 


Uhee vn = (AVin) + m, VmAo! Shale, 


rearranging terms yields: 
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The change in effective stress may be expressed as: 
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where: Ao = change in total stress 


The change in total stress at depth d may be represented 


as: 

Ag =. Ay;.:a 

OX, 

Do = (Vg/Vimn) bY, + (Vw/Vm) bY Vd + (V,/Vm) AY 


Since the unit weight of the sand grains does not vary 


with temperature or pressure (Ay, = 0) then: 


Ka = (Ve 7 VnlAy.d Seta / Ve) Arwa 


As a result, equation 3.19 becomes: 
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equation 3.19 can be rewritten as: 


Oe pees Oe ee ee or CAV Ree Vin (Vy OY O/ Vin +) VwAdw 


G/Vinli = My, Vpou-m.Vyou,- mov, Au Seat 


SC Can he Veen eV canon VAVn) a= wa var Plot VAT 


(section 3.1.1) equation 3.21 becomes: 


Ce Veo we emnOe A anv al  o Veal) “anime = My Ven Ve oye 


SEV TOT, Os= My Vga ome, awe me AU BA? 


Assuming the oil sand mass is fully saturated, the 
extent to which the voids are filled by bitumen or water 


(bitumen and water porosities) can be expressed as: 
n= V,/Vm and n= Vv/Vm 32123 
equation 3.22 may be rewritten as: 
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Ay, 0) / My t nem. th my Se 24 


It should be noted that a,, a,and ayare positive in 
the above equations since an increase in temperature 
results in an increase in volume; whereas, my, Mp and My 
are negative because volume decreases with increasing 


pressure. The" coefficient a,; is positive when an 
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increase in temperature produces an increase in volume 
of the sand structure and negative when a decrease in 


volume of the sand structure occurs. 


3.2 Analysis of Fluid Flow 


3.2. tintroduction 

To evaluate the rate at which excess pore fluid 
pressure resulting from an increase in temperature 
dissipates from an oil sand mass, the rate of fluid flow 
must be examined. One dimensional fluid flow through a 
Saturated soil 1S governed by Darcy's observational law 


which is represented by the equation: 


Mos +key Bai25 
or 

G = theraA 320 
where: 

v = flow velocity 


k = hydraulic conductivity 


he. 
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hydraulic gradient 


flow rate 


1Q 
uN 


A = cross-sectional area of soil normal to the direction of 
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flow 


The value of the hydraulic conductivity for a soil is 
influence by the properties of the permeant and the soil. 
The Kozeny-Carman equation provides a good assessment of the 
hydraulic conductivity in uniform sands (Mitchell, 1976) and 
illustrates the effect soil and permeant properties have on 
permeability. The Kozeny-Carman equation is expressed in the 


following manner 


Kae yerd/makos? abi tce) 852% 
where: 

k = the Darcy hydraulic conductivity 

Y = unit weight of permeant 

e = void ratio 


Ko = pore shape factor 


S = the specific surface per unit volume of particles 
mM = viscosity of permeant 

It 1S apparent from the Kozeny-Carman equation that the 
hydraulic conductivity is dependent, in part, on the 


viscosity (and Ptheeunit weight of ithe *permeant i9Temperature 
changes influence both the viscosity and the unit weight (or 
density) of the permeant and therefore the hydraulic 


conductivity also varies with temperature. It is useful to 
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express hydraulic conductivity in terms of absolute 
permeability K (a coefficient dependant on the properties of 


the soil skeleton only) by the equation: 
Rea Kk (yfAT)] 7 2fdty) o228 


Lambe and Whitman (1969) outline the soil 
characteristics which affect hydraulic conductivity. These 
characteristics include void ratio, particle size, fabric, 
composition and degree of saturation. When considering an 
oil sand mass subjected to a change in temperature changes 
in particle size, void ratio and fabric can be expected. 

Consider the time varient flow through an element of 
fully-saturated soil having dimensions dx, dy and dz with 
plow occurring inthe, ™,y oplane “only (Figure 3.1). The 


volume of fluid entering the element during a time At is: 
ive > ovy/ exidx/2)dyazaAret ily. = (dv, /oyiay/2)dxdzAt 
where: 


Vx, Vy = components of velocity at the central point P (x = 
y = z = 0) of the element 

Hevu/ox)), Love /dy) =echange of velocity in the « and -y 
directions 


dydz, dxdz = cross-sectional areas 
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Figure 3.1 Time Varient Flow Through an Element of Fully-Saturated Soil 
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The volume of fluid leaving the element during a time At is: 


Remedy dxsax7 2c ayazAt st (vy, + Lov, /oy Jdy/2)dxdzat 


Therefore, the net volume leaving the element during a time 


At ise 


Ohavy/ ox), + Lov, /ey)])dxdydzAt 


During the time increment At, the total head at point P 


increases by AH and as a result the volume of fluid taken 


into storage due to the change in total head is: 


S,dxdydzAH 


where: S, = specific storage coefficient 


The volume being generated within the element due to the 


thermal expansion of the pore fluid during the time At is: 


VgAt = g'dxdydzAt 


where: g' = thermal generation rate 


From the principle of continuity (conservation of volume) 


the above three quantities must sum to zero, hence: 
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(Lav,/dx] + [dv,/dy])dxdydzAt + S,dxdydzAH +g'dxdydzAt = 0 


substituting: 

v, = ~k,[dH/ax] 

and 

vy = ~k,[dH/dy] 

into equation 3.28 and rearranging terms yields: 
8/dx(k,[dH/dx]) + 3/dy(ky[dH/dy]) + g’ = S,[dH/at] 3.30 


Smith and Chapman (1983) identify two driving forces 
which causes fluid to flow in a thermal regime, piezometric 
head differences and Be ceyaacy force due to differences in 
fluid densities. When considering a heated foundation on oil 
sand, however, the direction of fluid flow is opposite to 
that of heat flow. As a result, the density of the pore 
fluid decreases ime ther tireccion (of!) eiluid):.flow mand 
therefore, the driving force associated with buoyancy 
effects does not exist. The effect of the fluid density 
decreasing towards the heat source as occurs in this problem 


1S included in the analytical procedure developed. 
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The generation rate (g') is the ratio of the amount of 
fluid expelled from an element of oil sand subjected to an 
increase in heat to the original volume of the element, 


pre .* 

la WV SAT! OeVyAT—ay AT (VT V La. eVmAT) /Vim REACH 

Since Vt Vy Hay =V,, 1, = Vy /Vm and ny= Vo /Vm 

aaa ™ oe AT se TagAle > a ATOny + Nw) — a -¢AT Beoe 
The specific storage coefficient, S,, is defined as the 

volume of fluid expelled per unit volume of aquifer in 

response to an unit decrease in total head. The pore volume 

compressibility, the compressibility of fluid, the density 

of fluid and porosity all contribute to specific storage. 

The specific storage coefficient can be expressed as: 


Se = CAV) (One Sei) CAV) PARP ST = (AV) (Ane ao119/V,) 93.33 


where: H = total head, which is the sum of the pressure head 


(u/y, ) and the elevation head (z) 


The change in total head may be written as: 
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AH = Au/y, + Az 


Since for a specific element of oil sand Az = 0 the change 


in total head is: 

AH = Au/y, 

For a unit change in total head (AH = 1) 

1-= Au/y, or Au = yz 

Since (AVw).5= m Vy Au 

then (AVy)[4H = 1] = myVyry 

and similarly 

(Av, ) [AH = 1] = mV, y,and (AV,) [4H = 1] = myVanrg 
As a result equation 3.32 can be rewritten as: 
Saoesdlm Vege + me Viige Be eA in 3.34 
Letting nj= V\/Vm and ny= W/Vm 


then, 
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Sig =engMyew PeNgMe %p ~My Ye SESE 

The above expressions allow the rate at which excess 
pore Fluid pressure resulting from an increase in 
temperature dissipates from an oil sand mass to be 
evaluated. When heat flow theory is coupled to equation 3.30 
(the fluid flow expression), the state of drainage or heat 


consolidation can be determined. 


3.2.2 Two Dimensional Transient Fluid Flow Beneath a Heated 

Foundation by Finite Elements 

The solution to the governing differential fluid flow 
equation for two dimensional transient fluid flow beneath a 
heated foundation (equation 3.30) was determined utilizing 
the finite element technique. Unlike the conventional 
governing differential equation representing transient 
seepage, equation 3.30 includes a term which accounts for 
generation of pore fluid due to increases in temperature. As 
a result, existing finite element programs designed to solve 
the conventional governing differential equation prove 
inadequate. Since heat transfer governing equations are 
compatible to the analysis of seepage problems and an 
anology can be made between the thermal generation rate 
(developed in this study) and the internal heat generation 
rate included in the heat transfer governing equation, an 
existing computer program which evaluates problems of heat 
conduction by the finite element method can be employed. 


When utilizing a heat transfer computer program to solve a 
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fluid flow problem, the appropriate analogous variable must 


be used. The analogy between variables include: 


temperature » total head 
thermal conductivity > hydraulic conductivity (permeability) 


specific heat » specific storage 


Using the above analogous variables, the heat transfer 
program ADINAT (A Finite Element Program for Automatic 
Dynamic Incremental Nonlinear Analysis of Temperatures) was 
employed directly for the analysis of the fluid flow problem 
(Chapter 6). Bathe (1977) provides a summary of the finite 
element theory as well as the ADINAT program. solution 
capabilities and the numerical techniques employed. 

This study determines the solution for fluid flow 
beneath a heated storage tank. Due to symmetry only one half 
of the tank is considered with no flow occurring across’ the 
axis of symmetry. The phreatic surface is assumed to be at 
the ground surface, below which the pore fluid is assumed to 
initially be static (hydrostatic conditions). The foundation 


material is also assumed to be fully saturated. 
3.3 Analysis of Heat Flow 
Sesedtdntroduction 


The conduction of heat from one mass to another or from 


one portion of a mass. to another is governed by two 
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fundamental laws of thermodynamics. The first law states 
that the thermal energy of a system is conserved (i.e., the 
thermal energy flowing into the system and the thermal 
energy generated inside the system must balance the energy 
flowing out of the system and _ the thermal energy stored 
inside th system). The equation which describes this concept 


is: 


Pam *Lq = Eout +E£s SIA Sis) 


where: 


Fin = rate of energy flow into the system 
Eg = rate of thermal energy generated inside the system 
Eout = rate of energy flow out of the system 


Es = rate of energy storage inside thesystem 


The second law of thermodynamics states that heat will 
flow only when a temperature differential exists and that 
Enem flow of heat, will takey (place from the point off the 
greatest temperature to the point of the least temperature. 
In other words, the flow of heat occurs only when a 
temperature gradient exists and occurs in the direction of 
diminishing temperature. The major processes of energy 
transfer are conduction, convection, thermal energy storage 


and energy generation. 
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through 


substances without motion of the conducting body as a whole. 


The rate at which heat flows by conduction can be expressed 


as follows: 


m= Kr AST /dx) 


where: 


k = thermal conductivity 


A = area normal to direction of heat flow 


1 = temperature 


x = distance 


Convection refers to the transference of heat between a 


solid and a circulating fluid. The rate at which heat 


by convection is given by the following expression: 


where: 


h = heat transfer coefficient 
T, = surficial temperature of solid 
Ts= ambient fluid temperature 


A, = surficial area 
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Thermal energy storage occurs when a solid is subjected 
to a time dependant increase in temperature. The rate of 


energy storage in a solid can be expressed as follows: 
ibee=lp.VLortatT/dt) 3239 
wheres: 


p = density 
V = volume 
c = specific heat 


@r=atime 


Thermal energy generation occurs in solids when other 
types of energy are transformed into thermal energy. The 
rate of thermal energy generation is commonly described by 


the following: 


where: 


thermal generation per unit volume 


< 
uN 


volume 


Performing an energy balance in two dimensions equating 


the time rate change of energy stored to the sum of net heat 
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transfer due to conduction and the heat generated yields the 


governing differential heat - conduction equation: 


cp(aT/dt) = 0/ox(k,,[dT/dx] ) + a/ax(ky, baT/ayd) + 
Poy Keyhorsexd dete o/ay (ky ys hdT/dyJddetog euchy 
The above equation is for non homogeneous’ and 


anisotropic heat flow, and simplifies for homogeneous, 
isotropic material where the thermal conductivity is 


eonstant,,t0o: 


Su /atwan(k7oed4 Lo? Share last hoaT/oyi dk #igzec 3.42 


The governing differential heat-conduction equation 
neglects the influence fluid flow has on the rate of heat 
flow. The amount of heat loss through convection is_ small 
considering that little fluid flows from the o0il sand 
Strata. Ignoring the effects of convection results in a 
Slight over estimation of the rate at which heat flows into 
the oil sand strata and the amount of pore fluid pressure 


and volume change that would occur. 


3.3.2 Two Dimensional Transient Heat Flow Beneath a Heated 
Foundation by Finite Differences 
The solution to the heat flow equation for two 
dimensional transient heat flow beneath a hot foundation was 


determined utilizing the finite difference technique. When 
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the second order derivatives of the heat equation are finite 
difference a system of ordinary differential equations 
(equations containing derivatives of a single independent 
variable) is obtained. A solution to this set of equations 
is determined numerically as a function of time. 

This study determines’ the solution for heat flow 
beneath a heated storage tank. Due to symmetry only one half 
of the tank is considered with no energy conducted across 
the axis of symmetry. The tank and ambient air temperatures 
are assumed constant at 175 and 10 °C, respectively. The 
initial foundation soil temperature was assumed to be 5 °C. 

The initial ee in the Finite difference formulation 
involves establishing a system of nodes which divide the 
area into a finite number of elements. A typical nodal 
arrangement for a portion of the region under consideration 
is illustrated in Figure 3.2. 

Figure 3.3 shows the energy system assigned to _ the 
interior node positioned at the coordinates (Xm, Ym) along 
with the appropriate energy terms. For the prescribed system 
energy flow occurs by conduction through each side and 
thermal energy storage within the system. 
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Figure 3.3 Energy System Assigned to an Interior Node 
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Figure 3.5 Energy System Assigned to a Node at Ground Surface 
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The rate equations may be approximately stated as: 
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The storage term may be written as; 
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Substituting the above rate equations into the energy 
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Nodes along the axis of Symmetry must be treated 
separately since no energy is transferred across the face 


(Figure 3.4), the energy balance therefore is: 


Geet. 6 Gort Gaiveet  * ESsieh 


The equations describing rate of heat flow are: 
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Nodes along the base of the tank and along the _ ground 
surface are also treated separately from the interior nodes 
due to convection across the boundary (Figure 3.5), the 


energy balance is therefore: 
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The storage term may be written as: 
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initial condttionk 


3.4 Coupling of Heat and Fluid Flow 

The expressions presented in section 3.1.1 allow the 
performance of an oil sand mass when heated to be evaluated 
under fully drained and undrained conditions. When heat flow 
theory (section 3.3) is coupled to the fluid flow expression 
(section 3.2), the state of drainage or heat consolidation 
can be determined; that is, whether the oil sand mass will 
be fully drained, fully undrained, or partially drained. If 
partial drainage occurs, the volume change and pore pressure 
change can be calculated. 

In this study a simplified approach was taken in the 
coupling of heat and fluid flow. The procedure adopted 
involves evaluating for each time increment the rate at 
which the heat from the structure propagates into the oil 
Sand strata, calculating the change in pore fluid pressure 
and change in engineering properties of the oil sand strata 
in response to the change in temperature, and subsequently 
determining the state of drainage in the oil sand strata. 
Complete coupling was not achieved since the influence fluid 
flow has on the rate of heat flow was neglected. The effects 
of convection, however, are small considering that little 
fluid flows from the oil sand strata. The coupling method is 


discussed in detail in Chapter 6. 
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3.5 Conclusions 

An expression was developed to calculate the net volume 
change of an oil sand mass subjected to a change in 
temperature under fully drained conditions. The volume 
change of the oil sand mass arises from the thermal 
expansion of the individual sand grains and structural 
changes in the sand skeleton. 

An equation which allows the volume change of oil sand 
subjected to increases in temperatures under completely 
undrained conditions to be calculated was presented. The 
volume change of the oil sand mass is due to expansion of 
the individual sand grains and pore fluids as well as 
structural changes in the sand skeleton. 

An expression was derived which expresses the volume of 
pore fluid which would drain from a fully saturated 011 sand 
mass experiencing a change in temperature. The amount of 
fluid expelled from an oil sand mass during a drained 
thermal expansion test is the sum of the volume change of 
the pore fluids less the change in void volume and 
Structural volume changes. 

An expression was developed to calculate the increase 
in pore fluid pressure in a fully saturated oil sand mass 
Subjected to a temperature change under fully undrained 
conditions. 

The theoretical expressions to evaluate the rate of 
fluid from heating were derived and are presented. The 


- method utilized to solve the governing differential fluid 
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flow equation for two dimensional transient fluid flow 
beneath a heated foundation using the finite element 
technique is described. 

The theoretical formulations to determine the rate of 
heat flow were derived and are presented. The solution to 
the heat equation for two dimensional transient heat flow 
beneath a heated foundation utilizing the finite difference 
technique is illustrated. 

The expressions developed in this study provide a 
method to analyze the rate at which heat flows through an 
oil sands mass, the rate at which excess pore fluid pressure 
resulting from an increase in temperature dissipates, and 


the amount of volume change which occurs. 
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4, LABORATORY APPARATUS AND PROCEDURES 


4.1 Introduction 

Inherent in the processing of oil sands at surface 
mining projects is the requirement of storage of heated 
bitumen in large retaining tanks. To experimentally 
investigate the problems which arise when these structures 
and other structures generating heat (e.g., hot water 
extraction plant, cokers, etc.) must be situated above oil 
sand strata, equipment has been built at the University of 
Alberta, Department of Civil Engineering. 

The. apparatus has the capability to. explore the 
one-dimensional thermal expansion of oil sands and its 
constituent components at temperatures approaching 200 °C 
and pressures in excess of 3.5 Mpa. One-dimensional testing 
(zero lateral strain, i.e., volume change due to expansion 
in the vertical direction) is justified considering that the 
minor principal stress in the surface mineable area is 
believed to be vertical from the surface to a depth of 
approximately 330 m and is one-third that of the horizontal 
stress to a depth of about 100 m (Dusseault, 1977b). 

Specifications and detailed dimensions of each 


component of the apparatus is given in Appendix B. 
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4.2 Laboratory Apparatus 
The following is intended to be a brief summary of each 


constituent of the equipment. 


4.2.1 Oedometer 

The assembled apparatus is expressed schematically in 
Figure 4.1 and is shown in Figure 4.2. The oedometer is 
conceptually simple in design. The sample jacket which 
provides lateral confinement for the specimen is held 
Stationary by a covering plate which is bolted to the base 
of the assembly. The specimen chamber is formed when the 
piston is lowered through the covering plate and is allowed 
to rest on the sample. This configuration ensures that 
vertical expansion or contraction of the specimen can be 
measured directly by recording the relative movement of the 
piston. Sealing of the sample chamber is accomplished with 
the aid of heat and petroleum resistant "O" rings, two of 
which are located on the lower part of the piston and one 
between the sample jacket and assembly base. A teflon rider 
ring is situated on the upper portion of the piston (Figure 
4.1) to avert tilting and binding during vertical movement. 

The oedometer has four drainage ports, two of which are 
contained within the piston and provide a drainage path from 
the top porous stone to the top of the piston, the other two 
provide a drainage path from the bottom porous stone and 
emerge through the base of the assembly below the sample 


chamber. One of the latter ports provides communication to 
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Figure 4.1 Drawing of Laboratory Apparatus 
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Figure 4.2 Laboratory Apparatus 
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the back pressure system. Located between the cell and the 
back pressure system is a stainless steel bottle, the 
purpose of which is to inhibit the flow of bitumen into the 
back pressure system. Valves are situated above and below 


the bottle. 


4.2.2 Axial Loading and Back Pressure Systems 

A diaphram air cylinder is used to apply the axial load 
to the oedometer piston. Briefly, air under pressure is 
forced into the cylinder causing the diaphram to expand 
against an internal piston producing the axial force. This 
force 1S transmitted to the oedometer piston through a small 
concrete cylinder, which retards the flow of heat from the 
oedometer to the air cylinder. 

The back pressure system (Figure 4.3) performs on the 
Same principal as the axial loading system. In this case, 
however, the diaphram in the accumulator 1S nitrogen driven 
againstesilicon ofl sutThe silicon oll int’ turn forces water 


under pressure to the oedometer. 


4.2.3 Heating System 

The heating system employed consists of three parts: a 
sensing thermocouple, a double coil electrical element band 
heater, anda digital temperature controller unit. The band — 
heater is clamped to the sample jacket as shown in Figure 


4.1. 
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The temperature controller unit (Figure 4.3) regulates 
the heat supplied to the sample chamber. The temperature 
controller once engaged will sense and display the 
temperature of the thermocouple. When a predetermined set 
point temperature is entered into the unit it will initially 
calculate and display the difference between the set point 
temperature and the temperature the thermocouple is sensing. 
If the set point temperature exceeds the thermocouple 
temperature the temperature controller will supply current 
to the elements of the band heater. Initially the cycle of 
heating iS rapid but as the set point temperature is 
approached the cycle time will reduce and when it:‘is reached 
the cycle time iS proportioned to hold the temperature at 
the thermocouple constant. In operational use it was 
observed that the time to reach the set point temperature 
was extremely rapid and stabilization occurred generally 


after only a few minutes. 


4.2.4 Measuring Devices - Data Logger System 

The measuring devices associated with the oedometer 
consisted of two pressure transducers, one linearly varying 
displacement transducer (LVDT ) and three Type J 
thermocouples. The specifications of each measuring device 
are discussed in Appendix B. 

One pressure transducer is situated to measure the air 
pressure in the diaphram of the air cylinder while the other 


is located in the line between the back pressure system and 


ee a 
Ae 


nagaluges 4 1 s ouein 
auld putas: adh, Cd | 
it elbegs a ean Seinen ‘Hie oe ’ 
yen Renin aioe a va | eigen Os: 


* rd f a 
VADBE D4 W 
' 4 
* , fo ea 
io Sey any = 


> 
< 
co) 
ty 
% 
rary 
= 


a. sae Saaes peed st $ verges tegmie st Wii 


i et ae shies 748) Pena sreqiioe “gata | 


F : ‘ : j ri e ‘ PASS : | / / | 
baa Views 6liw’ Pers a ane ymesmbd ne 


ci Syadayediio. faier ri Be soak Sige. C7 
baldipaw? at (4 soiw hag obey ate ant ase: 
He Biv eoehes, oft Blok, od” ‘ote narsoerone at mh, 
x ¥ i whe Lane eS gi nit er WO, eh 


Tusa rAmes hae awe oni Acme, a ai aa sed 


AS 
i lgae tee: et wsd¢ | bet hres es bee. biyer el 


fe estab et & r. 
PNR ; 
‘dean evap eed <i aes tsi aby 


Wy 
4 


‘dihah Se eRe (hos i aed ao8zs | an ivie6 ‘pilsoehun 


nerieaey vi" RSI S04; etaoubenets ayyees1g .ovs 
be peer nea (Tava) - “Ragubans ya 


« i oe 


O51¥ee. Poluvecam  ‘doay to ‘Sao nt sapisgede vet 


2 ePonsnyA' 2 
ats. abe Sunken -o- nee ai i Nei ee | 
wadto ii S4 dw robab ys Tis ons tocis mabey ts UB eit? ne 


besa ie2 Siveneny ae ore rer PEL) acts ar: 


iad fig ; fr 
3 aa ; va ; ae *} 7‘ Pas 
= 7 : Pts a 
oh | me th eee ies eae yy mike 
‘ , la! Ye : 
i] p Ve q 7 
; ve ah 


61 


the oedometer. One end of the LVDT is secured to the frame 
of the apparatus below the oedometer and the other end to 
the piston rod of the diaphram air cylinder. The 
thermocouples are Situated as follows: one in the oedometer 
piston with the meaSuring tip in contact with the top porous 
stone; one at mid-sample height in the sample jacket; and 
one in the base of the apparatus in contact with the bottom 
porous stone (Figure 4.1). 

The data logging was performed by a Hewlett - Packard 
3497A Data Acquisition Control] Unit which electronically 
gathered and recorded the data during the testing program 


(Figure s4.4).. 


4.2.5 Insulation 

The base of the oedometer iS mounted on sheets of 
asbestos to reduce heat loss to the steel frame of the 
apparatus (Figure. 4.j)agourang@soperationad, use , the 
oedometer itself is completely enclosed in an insulating 


shell composed of ceramic fiber insulation (Figure 4.5). 
4.3 Test Procedures 


4.3.1 Thermal Expansion of Water 

Two thermal expansion tests were carried out on 
distilled water under pressures of 2000 and 3500 kPa. In 
both cases the sample chamber was filled at room temperature 


with water to a height of approximately 2.5 cm. All ports, 


wmnat ads or 


Deen” DSA ae Ma 
; 


t 
: Endl 
i: Tae. TA 
ay 
2 LOL J 
ne Se + 


OF iS 
WW end 4 


near 


maid 


tie 


cwlaine Riso has +a Bon ont 
4 Be 


; ee 


Boo oes vollied 


wag 


venta wiobeaaisip bee 


ht aon) Wel 2 ae aap boa atgnee “ot ee 
i ped ais ew ng aig “nit ue ry ayes “ota ie iiss 
a —7 1 , | 
5 m ‘: 7 a) 
i it t #4 wa 
3 sete Wo Baeomia a! ot wigp bed ‘ata a 
Wo pam yee: SN) Ru eat sae at 
joy, TahoLje age grkaue, oes rt 
: 1 sv auitceties ae a sha 


pa 


vor bay ehy 9e 


eo a 
“pt wm | elnsinth 
ee sedate Sigass sts cone dad” ‘i 


[a alain, 


sateen) Gebhe 


shan! vonage tee dipian s o2, vosaw' date! 
ss Pir , te 


62 


Figure 4.4 Data Logger System 


Figure 4.5 Laboratory Apparatus Enclosed 


in Insulating 


Shell 
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passages and porous stones were saturated prior to testing. 
Each test was performed under a constant back pressure. 
The back pressure was held constant during the test through 
manipulation of the applied load to the piston. The heating 
increment usediVetor. Gach. itest: \ewasi) 25)" “.@yeeAtter sceach 
application of heat, sufficient time was permitted to 
establish steady state conditions, i.e., the temperature 
readings had stabilized and the thermal expansion (volume 


change) ceased. 


4.3.2 Thermal Expansion of Bitumen 

Seven thermal expansion tests were performed on bitumen 
under pressures of 500 to 3500 kPa and temperature range of 
22 to 200 °C. In all cases the sample chamber was filled at 
room temperature to a height of approximately Zo 
centimeters with bitumen and all ports, passages and porous 
stones saturated with distilled water. Prior to each test 
the system was permitted to saturate (gases driven into 
solution) for a period of up to 24 hours at the appropriate 
pressure level. 

Each test was performed under a constant back pressure 
with heating increments varying from 5 to 25 °C. After each 
application of heat sufficient time was permitted to 
establish steady state conditions, i.e., the temperature 
readings had stabilized and the thermal expansion (volume 


change) ceased. 
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4.3.3 Thermal Expansion of Compacted Tailing Sand 

Four tests were performed to investigate the thermal 
expansion of tailing sand at various densities. The 
temperature range for the testing was varied from room 
temperature to 200 °C. To eliminate any influence of induced 
pressure within the sample chamber, the sand was tested in 
Piesgatr dried.pstate,eand «during. ~thex, testing) the, bottom 
drainage ports were held open to atmospheric pressure. 

The initial test investigated the thermal expansion of 
sand in the loose state. To achieve a loose sample, sand was 
Slowly poured into the sample chamber and carefully levelled 
uSing a small metal spatual. The piston was then | placed 
inside the sample jacket and allowed to fall under its own 
weight to rest on the sample. Care was taken to avoid 
jarring of the apparatus. In order that the LVDT could be 
utilized to measure the vertical expansion (or contraction) 
of the specimen during the test, it was necessary to apply a 
small axial seating load of 50 kPa. The density obtained was 
1.46 Mg/m?. Two subsequent tests were performed on sand 
under the same seating load of 50 kPa but at higher 
densities. Two procedures were utilized in achieving a dense 
Sample. The first involved densifying the specimen by 
cycling the applied load from 50 to 3500 kPa. The second 
method involved placing the sand in the sample chamber in 
four lifts, vibrating between lifts under an applied load of 
50 kPa. The methods produced densities of 1.57 and 1.62 


Mg/m? respectively. 
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The pinwiadnernar expansion test was performed on sand 
under an applied load of 3500 kPa. The sand was densified in 
the sample chamber by vibrating in four lifts. The density 
obtained was 1.65 Mg/m°. 

The test procedure employed was the same for all four 
samples. It involved heating the sample in 25 °C increments 
allowing sufficient time for temperature stabilization to 


occur at each desired level. 


4.3.4 Thermal Expansion of Oil Sand 


4.3.4.1 Sample Preparation 

The oil sand specimens used in the testing program 
were obtained from 94.5 mm core samples taken from an 
oil sand outcrop st Saline Creek, approximately 1 km 
south of Fort McMurray, Alberta. The oil sand core 
Samples prior to sample preparation were sealed in 
plasticpand stored in-ascold room at -—20"ce +25: °C. 

When an oil sand sample was required for testing a 
piece of core was taken and sealed in an air-tight metal 
container. The container was then placed in a styrofoam 
receptacle and covered with dry ice. Sealing the sample 
in the metal ‘container prevents disturbance of the 
Sample from the carbon dioxide gas emitted from the dry 
ice. The container was left in thewedry ice 
(approximately -78 °C) for several hours. 

To? obtain: «a test Sample of the proper dimensions 


(7.62 cmmin® Gbameters approximately 2.5 cm in height) 
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the following procedure was employed (between each step, 
to inhibit swelling of the specimen from heat developed 
during machining, the sample was sealed in the metal 
container and placed in dry ice for a minimum of 30 


minutes): 


e The piece of core was placed in the lathe (situated in 
thescold roem)), “gumppedsat ionemendaby themmjawc, of the 
chuck and supported at the other end by a freely 
rotating disc (tail stock). The surface of the specimen 
at the end supported by the disc was machined using a 
tungesten carbide bit until round and smooth. The 
Specimen waS removed from the lathe and an adjustable 
metal clamp (2.5 cm wide) placed on the trimmed section. 
e The specimen was placed in the lathe with the end 
having the metal clamp gripped in the chuck. Reference 
marks were placed on the chuck and metal band. Similarly 
at the opposite end a reference mark on the support disc 
was aligned with a notch on the specimen. The surface of 
the sample was then machined from the metal band “to 
within 2.5 cm of the opposite end (Figure 4.6). Ten 
passes were made, reducing the diameter by approximately 
0.5 mm each time. 

e The sample was returned to the lathe, realigned with 
respect to the reference marks and 10 more passes made. 
e The previous step .waS repeated until the correct 


diameter was reached. 
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e The specimen was returned to the lathe, aligned and 
Square notches (about 2 cm deep) cut at 2.5 cm spacings 
(corresponding to the sample height). The notching 
ensured the sample ends would be parallel and square. A 
notched specimen is shown in Figure 4.7. 

e The Specimen was cut with a diamond saw at the 
notches. 

e The test sample supported by a split metal ring was 
placed in the ere and the central portion (cut by the 


saw) at each end machined flush (Figure 4.8). 


4.3.4.2 Undrained Test Procedure 

The prepared oil sand sample was inserted into the 
sample jacket in the cold room. The sample jacket was 
then transferred to the laboratory where the apparatus 
was quickly assembled and a 100 kPa seating load placed 
on the sample. The sample was allowed to thaw for 
approximately 4 hours after which time water was allowed 
to percolate up through the sample to saturate the 
passages and ports in the piston. The piston pressure 
was increased to slightly above the required test 
pressure. The back pressure was then increased 
incrementally to the desired level. The system was then 
allowed to saturate for about 24 hours to ensure all 
gases were dissolved in the liquid phases. 

The test was performed under a constant back 


pressure. The back pressure was held constant during the 
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test through manipulation of the applied load to the 
piston. The heating increment used for the test was 20 
*CasAEter each.application of, heat, sufficient, time, was 
permitted to establish steady state conditions, i.e., 
the temperature readings had stabilized and the thermal 


expansion (volume change) ceased. 


4.3.4.3 Drained Test Procedure 

The methods employed to mount and saturate the test 
Sample were the same as those adopted for the undrained 
test (Section 4.3.4.2). 

The drained test was performed under a. constant 
effective confining pressure (difference between back 
pressure and piston pressure) of 240 kPa. A volume 
change measuring device was inserted between the back 
pressure system and the fluid trap to electronically 
monitor the fluid drainage during the test. The heating 
increment used for the test was 20 °C. After each 
application of heat, sufficient time was permitted to 
establish steady state conditions, i.e., the temperature 
readings had stabilized and the thermal expansion 


(volume change) ceased. 


4.4 Conclusions | 

Equipment was successfully designed and constructed 
with the capability to explore the one-dimensional thermal 
expansion of oil sands ie its constituent components to 


temperatures approaching 200 °C and pressures in excess of 
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3.5 MPa. A brief summary of each constituent of the 
equipment is presented. 

The test procedures which were developed for use in the 
laboratory program worked successfully and are described. 
The method utilized to prepare oil sand samples from frozen 
core resulted in very high quality test specimens. The 


procedure adopted is illustrated. 
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5. EXPERIMENTAL RESULTS 
5.1 Thermal Expansion of Water 


5.1.1 Introduction 

Due to the amount of emphasis placed on data obtained 
in the testing program, it is critical to ensure that the 
experimental information derived from the apparatus proves 
precise and reliable. To establish this confidence in the 
apparatus, it waS necesSary to test a substance whose 
behavior under temperature and pressure is well documented. 
Water was selected as the testing material owing to the vast 


amount. of information contained in the Literature. 


5.1.2 Test Results 

Two thermal expansion testS were performed on water 
under pressures of 2000 AAG 3500 kPa and temperature range 
pees to 200 °C (Figures 5.1 and 5.2). For completeness, sam 
data gathered from the tests are plotted on the figures. It 
should be noted that steady state conditions for a 
Particular temperature level were only reached when the 
temperature readings had stabilized and the thermal 
expansion (volume change) ceased. This condition is 
reflected on the plot where there is a high density of data 
points. 

For comparative purposes theoretical curves were 


calculated based on information obtained from the 1968 
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edition of the vapanese Society of Mechanical Engineers 
Steam Tables. These curves reflect steady state conditions 
and are drawn on the plots. 

Based on the steady state criteria (stabilized readings 
for each temperature increment) there is excellent agreement 
between theoretical and experimental results in both cases, 
indicating the apparatus performance and calibrations are 
Satisfactory. 

The data reduction calculations, data reduction 
computer programs and experimental results are contained in 


Appendix D. 


5.2 Thermal Expansion of Bitumen 


5.2.1 Introduction 

The oil sands present in the mineable portions of the 
Athabasca Deposit have bitumen saturations ranging from 0 to 
18 weight percent (up to approximately 36 volume percent 
[Mossop, 1978]). Present within this bitumen- may be 
appreciable quantities of free and dissolved gases (mainly 
methane and carbon dioxide [Dusseault and Morgenstern, 
1977]). The hot water extraction method is utilized to 
Separate the bitumen from the other oil sand constituents. 
Briefly, the process involves subjecting oil sand to hot 
water, steam and caustic and. collecting the bitumen as 
£roth: Subsequent centifuging completes the separation 


process. The bitumen obtained from this process can be 
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expected to be relatively free of dissolved gases mand to 
have a slightly higher density (resulting from the lighter 


components being driven off during the process). 


5.2.2 Test Material 

Since no technology exists which allows undisturbed 
sampling of pure bitumen, bitumen obtained from the Syncrude 
extraction facility was used as the test material. Prior to 
testing, the bitumen was stored at room temperature ina 


sealed container. 


5.2.3 Test Results 

The data reduction Calculations, data reduction 
computer programs and experimental results are contained in 
Appendix D. Points representing steady state conditions for 
each pressure and temperature level are shown in Figure 5.3. 
Figure 5.4 provides a comparative plot of the experimental 
thermal expansion of bitumen and the theoretical curves for 
water under pressures of 500 kPa and 3500 kPa (as given in 
the previous section). Figure 5.5 gives a comparative plot 
of the density change of bitumen with temperature with that 
of water. 

Figure 5.3 indicates that the cumulative volume change 
of bitumen over the temperature range of room temperature to 
200 °C is not dependent on pressure. The thermal expansion 
of the bitumen tested is greater than that of water up to a 
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wcmperacurer range oLrr25rto"200=°cm(PigurerSedeqm gure 5.5 
Shows that the density change of bitumen with temperature is 
Similar to water from room temperature to about 125 °C and 
decreases at a slightly slower rate from 125 ° to 200 °. 
Figure 5.5 also shows that the bitumen had the same initial 


density for all pressure levels considered. 
5.3 Thermal Expansion of Compacted Tailing Sand 


o-56| Introduction 

Laboratory tests were performed on compacted tarling 
sand to investigate the overall thermal expansion of the oil 
sand skeletal grains. The measured thermal expansion is the 
combined volume change of the soil Structure and the 


individual sand grains. 


5.3.2 Test Material 

The sand used in the testing program was obtained from 
the Suncor operation and consisted of light brown, uniform, 
Sub-angular quartz sand. The grain size analysis produced a 
distribution shown in Figure 5.6. Following the recommended 
A.S.T.M. Procedures (D 2049-69) the maximum and minimum 
densities were found to be 1.65 Mg/m* and 1.41 Mg/m°, 
respectively. The minimum void ratio (@min waS 0.61 and the 
maximum void ratio (emx) 0.87. 

Hardy (1974) reports testing on sand tailings yielding 


a void ratio of 0.60 for sand in its most dense state and a 
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Figure 5.6 Grain Size Distribution of Oil Sand Tailings 
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void ratio of 0.87 in the most loose state. Hardy (1974) 
also provides two curves which encompass the majority of 
grain size distributions found in laboratory tests. For 
comparative purposes, these curves have been drawn on Figure 
5.6. Hardy (1974) describes the sand grains as sharp and 
abrasive. 

The visual properties and densities are consistant with 
those reported by Hardy (1974). Comparison between grain 
Size distributions indicate the sand used in the testing 


program contained less silt size material. 


5.3.3 Test Results 

The data reduction Calcularions: data reduction 
computer program and experimental results are contained in 
Appendix D. Points representing steady state cumulative 
Petane change at each temperature level (with respect to 
room temperature) for each particular density are shown in 
mequne- 5,7." Figure, 95.7 also ‘contains; for compare. 
purposes, the cumulative volume change curve for an alpha 
Quartz ecrystals A typical sand grain may be composedvofta 
number of randomly orientated quartz crystals and would 
expand ethermally with ‘no voids; forming. Thenetore; che 
thermal volume change of a sand grain can be approximated by 
the expansion of a single quartz crystal. Fagure. 5.8 
contains plots of density change versus temperature for each 


Sample tested. 
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For each thermal expansion test conducted, a negative 
volume change was observed for the first heating increment 
(room temperature to 50 °C [Figure 5.7]). The Joose sand 
sample (initial density of 1.46 Mg/m* at room temperature) 
displayed the greatest volume decrease (0.36 %) over this 
temperature range. The sample initially compacted toa 
density of 1.62 Mg/m? at room temperature showed a volume 
decrease of 0.08 % over the temperature increase from room 
temperature to 50 °C. After completion of this test (to 
approximately 200 °C) the sample was allowed to cool. At 
room temperature a volume decrease of 0.10 % was measured 
relative to the original volume at room temperature. 

Since the sand was tested in the dry state with the 
grains Lmawgnineral to ‘mineral contact and “each (grain 
undergoing thermal expansion, relative movement Or 
reorientation of sand grains must have occurred during the 
temperature increase to 50 °C to allow an overall volume 
decrease. Campanella and Mitchell (1968) suggest such 
behavior results from a temperature induced change in 
interpartical forces that occurs to allow the sand structure 
to support the same effective stress. 

Over thesitemperatune parangeiOf, SQer2Geto 200 —%Gsad) 
Samples display a volume iperease roughly paralleling that 
Geka quartzea crystal... (Rigures5.47) waThis dbehavionstends ato 
Suggest theysoilsistructure - has sachievedytae tight packing 
arrangement (grain slippage no longer ,,oeccurning sor 


negligible) and the volume increase with temperature result 
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from thermal expansion of the individual grains. 

Figure 5.8 illustrates the change in density of each 
Sand specimen with temperature. All samples tested display 
the same general trend, the density increaseing over the 
temperature range of room temperature to 50 °C and gradually 
falling off to below the original density (at room 
temperature) over the temperature range of 50 °C to 200 °C. 
The sample originally at a room temperature density of 1.62 
Mg/m? shows a slight increase in density upon cooling back 


to room temperature. 
5.4 Thermal Expansion of Oil Sand 


5.4.1 Introduction 

Oil sand subjected to increases in temperature will 
experience volume increases, the magnitude of which is 
dependent on the amount of pore fluid drainage which is 
permitted. The maximum thermal expansion occurs’ under 
undrained conditions and the minimum under conditions of 
full drainage of pore fluids. 

Thermal expansion tests performed on oil sand _ under 
undrained conditions generates volume increases of mineral 
grains, pore fluids and mineral skeleton. Thermal expansion 
tests carried out on oil sand under drained conditions 
causes volume increase of mineral grains and mineral 
skeleton. In either case, if adequate back pressure is not 


provided, exsolution of gases dissolved in the liquid phases 
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Willie. OCCULT 


5.4.2 Test Material 

The oil sand samples used in the experimental testing 
were obtained from core samples taken from an oil sand 
outcrop at Saline Creek, approximately 1 km south of Fort 
McMurray. The pertinent test sample data is given in table 


Dye al cs 


5.4.3 Test Results 

The results of undrained and drained thermal expansion 
tests performed on av igsata under a constant back pressure 
of 2000 kPa (confining pressure of 2240 kPa) are given in 
Figure 5.9. Also contained within this figure is a plot of 
volume of fluid drained from the sample and a curve of the 
combined volume of fluid drained and volume change of the 
drained sample. Figure 5.10 illustrates the change in bulk 
density with temperature for both the undrained and drained 
cases. 

Two additional drained thermal expansion tests were 
performed on Saline Creek oil sand under a back pressure of 
5000 kPa (confining pressures of 11000 kPa and 5200 kPa). 
The results, of. these tests: are given in Figures) S211. Also 
Contaaned. in thst figures asothe drained3test<carpiedvout 
under a back pressure of 2000 kPa (confining pressure of 
PatQeskPaye Eigure 5.12 gives the results of the drained 


thermal, expansion tests» in terms of dry density. For 
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Table 5.1 Test and Sample Data 


Type of 
Zest 


Back 
Pressure 
(kPa) 


Vertical 
Confining 
Pressure 
(kPa) 


Effective 
Confining 
Pressure 
(kPa) 


Bulk 
Density 
(Mg/m° ) 


Percent 
Bitumen 
by Weight 


Percent 
Sand 
by Weight 


Percent 
Water 
by Weight 


Undrained 


2000 


2000 


22027 


84.93 


Test 


Drained 


2000 


2240 


240 


Zea cel 


84.93 


Void Ratio 0.539 0.539 


Dry Density 1a fe few el 
(Mg/m? ) 


Number 


Drained 


5000 


11000 


6000 


WeSies 


Was 25) 


82.34 


89 


Drained 


5000 


5200 


200 


12.0 
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comparitive purposes, the results of the tests performed on 
compacted tailing sand (Section 5.3.3) are included in the 
figure. 

The drained thermal expansion test carried out on oil 
sand under a back pressure of 2000 kPa (Figures 5.9 and 
Bend ite displays a rapid increase in volume over the 
eenperagurerrangesofie 180r° Cato 200° Gea Priortko: 180an’Gerthe 
rate of volume increase with temperature is consistant with 
the results of the two other drained tests (Figure 5.11). 
This observation tends to suggest the back pressure was not 
sufficiently high enough over this temperature range to keep 
gases in solution. If gas exsolution did occur over this 
temperature range, the gas bubbles at time of formation must 
have been at pressures exceeding the vertical confining 
stress (2240 kPa) in order to disrupt the skeletal grains 
and force the piston upward. The gas bubbles would 
theoretically continue to expand until pressure equilization 
between the fluid and gas phases was achieved resulting in 
the steady state points (temperature stable and zero volume 
change) at 180 and 200 °C. 

The undrained and drained thermal expansion tests 
(under back pressure of 2000 kPa) were performed on similar 
oil sand samples obtained from the same core specimen. This 
implies the measured volume change of the drained sample 
combined with the volume of fluid drained from the sample 
should theoretically approximate the results of the 


undrained test (Section 5.4.1). The curve reflecting the 
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combined volume of fluid drained and volume change of the 
drained sample, although approximately parallel to the 
undrained test results, is consistently lower. The 
difference is due primarily to the contraction of the pore 
fluid as it leaves the sample chamber and is collected and 
cooled at room temperature in the fluid trap. The volume 
change apparatus does not take into account the contraction 
of expelled fluid and as ae result consistently under 
estimates the amount of fluid drained during the test. The 
undrained test in contrast to the drained test shows no 
evidence of gas exsolution. 

Several observations are made from the drained thermal 
expansion tests performed on oil sand. For the tests 
conducted under vertical confining pressures of 5200 and 
11000 kPa negative volume changes were observed after the 
initial heating increments (Figure 5.11). Similar behavior 
was observed for the tests on compacted sand (Figure 5.7). 
Figure 5.12 shows the dry density change with temperature of 
these drained tests are consistant with those of the 
compacted sand. The initial densities were also similar 
indicating appreciable swelling occurred during sample 
preparation. The drained test performed under 2000 kPa back 
pressure did not display contraction during application of 
initial heat increments and had a_ significantly higher 
ifitiahaduy density we(1s72reMg/mi) ceflectingr less” tsample 
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Also shown on Figure 5.11 is the cumulative volume 
change with temperature of the quartz sand grains. Since the 
result of a fully drained thermal expansion test on oil sand 
is a measure of the volume change of both the sand grains 
and soil structure, the difference between the two curves is 
the volume change of the soil structure (Figure 5.11). A 
plot of the cumulative volume change of the soil structure 
with temperature for the drained tests performed is given in 
Figure 5.13. 

An assumed cumulative volume change with temperature 
curve for in situ oil sand (average bulk density of 2.12 
Mg/m*; dry density of 1.80 Mg/m? [Dusseault, 1977]) is given 
Poe er igure: 5.13) Ew is sevident from /Figure, Deals mthat 
temperature induced changes in the. soil structure is 
Significant over the temperature range of room temperature 
to 140 °C and is greatest at approximately 50 °C. The more 
disturbed samples (bulk density of 1.968 and 1.990) showed a 
higher degree of contraction of the soil structure. 

Figure 5.14 gives the change in the incremental 
coefficient of thermal expansion of soil structure with 
temperature. The incremental coefficient of thermal 
expansion for oil sand tested increases rapidly over the 
temperature range of room temperature to 50 °C and becomes 
Somewhat constant over the range of 50 to 200 °C. The amount 
of change in the incremental coefficient of thermal 
expansion of the soil structure during the initial heating 


increments (to 50 °C) is more pronounced for the disturbed 
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samples (bulk densities of 1.968 and 1.990 Mg/m?). 

From the drained. test conducted on oil sand it. is 
evident that volume increase of less than 1 percent is 
likely for drained conditions over the temperature range of 
20 to 200 °C provided adequate back pressures are maintained 
to prevent gas exsolution. The thermal volume increase 
appears to be independent of pressure (Figure 5.11). For 
Saline Creek oil sand significant thermal volume increases 
would be expected under back pressures below 2000 kPa. 

Based on the undrained test results thermal volume 
increases of approximately 6 percent would be expected when 
Saline Creek oil sand is heated to 200 °C (assuming the back 


pressure is sufficient to keep gases in solution). 


5.4.4 Theoretical Analysis to Predict the Undrained Behavior 
of Oil Sand 
A theoretical analysis for predicting volume changes of 
oil sand due to temperature variations was presented in 
Chapter 3. From the theoretical analysis the change in 
yolumes-.O£i.0i1l+»sandsiunder undrained conditionsa,.ands-at 
constant effective stress was determined to be: 


(AVin GeVaATotya Ve AT tgea4V. ATata(aV, ode 


es = 
The above formulation is used to approximate the results 
obtained from the undrained thermal expansion test carried 


out on Saline Creek oil sand. 
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The coefficients of thermal expansion for bitumen and 
water used in the analysis was obtained from tests performed 
at the same back pressure as the undrained oil sand test 
(2000 kPa). Data collected from the drained test performed 
on oil sand obtained from the eame core sample and tested at 
the same back pressure was used as input for the sand grain 
and sand structure terms. Figure 5.15 shows the experimental 
results of the undrained thermal expansion test along with 
the results of the theoretical analysis. 

The theoretical analysis provided an underestimation of 
the thermal volume change obtained from the undrained test. 
Since the results obtained from the thermal expansion tests 
performed on water displayed excellent agreement with 
published data and results of the drained test (with regard 
to behavior of the sand grains and sand structure) were 
consistant with thermal expansion tests on Similar oil sand 
(as well as tailing sand) then questions arise as to the 
validity of assuming results obtained from tests performed 
on extracted bitumen are representative of in situ bitumen. 
Too as feasible that bitumen extracted» ‘by -the» her water 
extraction process may have undergone some physical and 
chemical alteration. In addition, the compositional nature 
of bitumen: “mayo ovary with respeety to locaton Gaeq, 
depositional enviroment). 

it it is assumed that extracted bitumen is not 
representative of in situ. bitumen, then results obtained 


from undrained and drained tests can be used as input in the 
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theoretical analysis to generate a curve representative of 
the thermal expansion of in situ bitumen. Such a curve was 
developed for Saline Creek oil sand and is shown in Figure 
5.16. For comparative purposes the results of the test 
performed on the extracted bitumen is also shown in _ the 
figure. The results indicate that in situ gas saturated 
bitumen may show temperature volume increases 50% greater 


than the extracted gas free bitumen. 


5.5 Conclusions 

Thermal expansion tests performed on water resulted in 
volume increases of about 15 per cent at 200 °C. The thermal 
volume change of water did not vary significantly for the 
pressure levels employed. Excellent agreement between 
experimental and published data indicates the pressure cell 
employed was calibrated and operating properly. 

Thermal expansion experiments carried out on hot water 
extracted bitumen resulted in volume increases of about 15 
per cent at 200 °C. The test results indicate that the 
cummulative volume change of bitumen over the temperature 
range considered is not pressure dependent. The volume 
change of this dead bitumen is more linear with temperature 
than that of water. 

For each thermal expansion test conducted on compacted 
tailing sand, a negative volume change was observed for the 
first: heating increment droom*temperature to 508 °C). The 


loosest sand sample (initial density of 1.46 Mg/m*® at room 
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temperature) gave the greatest volume decrease (0.36%) over 
this temperature range. The sample initially compacted to a 
density of 1.62 Mg/m® at room temperature showed a volume 
decrease of 0.08% over the temperature range from room 
temperature to 50 °C. After heating to 200 °C, this sample 
was cooled to room temperature anda volume decrease of 
0.10% was measured indicating that the volume decrease was 
due to a compaction of the soil structure. Relative movement 
or reorientation of Sand grains must have occurred during 
the temperature increase to 50 °C to allow an overall volume 
Secreasce. Over =the temperature range ot 50 to 5200)" °Cs all 
Samples tested displayed a volume increase roughly paralling 
@oaemoOr a Quartz. crystal, indicating the ~so1l "structure 
achieved a tight packing arrangement at 50 °C and the 
Subsequent volume increase with temperature results only 
from the thermal expansion of the individual grains. 

Undrained heating of oil sand to 200 °C will result in 
a volume increase of about 6 percent. Most of this increase, 
over 5 percent, is due to the pore fluids, water and 
bitumen, with the bitumen expanding slightly more than 
water. The remainder of the volume increase is due to _ the 
thermal expansion of the sand grains. The dense, 
impenetrative structure»of in situ. o2l. Sand® “will show) a 
negligible amount of thermal volume change by itself. 

From the tests performed on oil sand the thermal volume 
increase appears to be independent of the effective 


confining pressure for the pressure range considered in this 
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Study. For the oil sand tested here, significant thermal 
volume increases would be expected under back pressures 
below 2000 kPa because of gas exsolution at temperatures 
above 180 °C. 

Results obtained from the undrained and drained thermal 
expansion tests performed on oil sand were used to generate 
a plot representative of the thermal expansion of in = situ 
bitumen. The results indicate that in situ gas saturated 
bitumen may show temperature volume increases 50% greater 


than the extracted gas free bitumen. 
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6. NUMERICAL ANALYSIS 


6.1 Introduction 

To determine the amount of heave which would occur 
beneath a structure where the oil sand foundation 
experiences a temperature increase, it is necessary to 
assess the spatial temperature increase with time and 
estimate the amount of drainage which might occur. Since 
such field observations are absent in the published 
literature, the rate of heat flow beneath the heated 
foundation and the amount of drainage which occurs are 
evaluated using numerical methods for a general 
Stratigraphy. 

The analytical model put forth here investigates the 
heat flow from an 80 m diameter storage tank resting on a 2 
m crushed rock and sand foundation pad. The tank is assumed 
to be held at a constant temperature of 175 °C, The general 
Stratigraphy chosen for one study consists of the granular 
Bad on top of «32 mot clayey till overtying rich oil sand: 

The analytical model developed in this study couples 
heat flow theory to the fluid flow expression to assess the 
State of drainage in the oil sand strata below the heated 
foundation and calculates the change in pore pressure and 
the volume change (or heave). The analytical model consists 
of four computer programs. FEHT2-1, developed by G. E. Myers 
atethe: University tot Colorado, was used tosevaluate the rate 


at which the heat from the foundation propagates through the 
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oil sand strata. VARTEST5, developed by the author, was used 
to calculate the change in pore pressure and change in 
engineering properties of the oil sand strata in response to 
the change in temperature (as predicted by FEHT2-1). 
VARTEST5 then compiled the information and formated the data 
Pots airect@ input’ into’ the-"fluid® flow programe? in*other 
words, VARTEST5 served to couple the heat flow analysis to 
the fluid flow analysis. ADINAT, developed by K. J. Bathe at 
M.I.T., was used to determine the state of drainage in the 
oil sand strata. The above three computer programs were 
excecuted in succession at predetermined. time intervals to 
Simulate simultaneous heating and draining of the oil sand 
Scratas 

A fourth computer program, VOLTEST5 (also developed by 
the author) was utilized to determine the amount of volume 
change occurring after each heating and draining cycle and 
to calculate the amount of heave occurring under the 


foundation. 


6.2 Heat Flow 

The thermal analysis was performed utilizing a computer 
program (FEHT2-1) developed by G. E. Myers at the University 
of Colorado in 1972 (Appendix E). FEHT2-7 (Finite Element 
Heat Transfer - 2 dimensions) evaluates problems of heat 
conduction in two dimensions by a finite element method. The 
method involves the use of linear triangular elements. The 


computer program has the capability of handling steady and 
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transient states. The transient states can be analyzed by 
any Of three methods: Eulers; Crank-Nicholson; or 
pure-implicit. The program considers three boundary 
conditions: time independent convection; time independent 
heat flux; and time independent temperature. The program was 
Originally dimensioned to accommodate 50 nodal points but 
has been altered to handle in excess of 280. 

The finite element mesh utilized for the numerical 
analyses is shown in Figure 6.1. The finite element mesh 
consists of 500 triangular elements and 284 nodes. The 
elements were proportioned in size to provide greater 
accuracy in critical areas. 

The Crank-Nicholson method (a central finite difference 
scheme) was utilized to determine the transient temperature 
variation beneath the Neneed foundation. Myers (1971) 
provides an in-depth discussion of the Crank-Nicholson 
scheme. The transient temperature variation was calculated 
using FEHT2-1 based on a time increment of 12 hours for the 
first Year and a time increment Of G0 hours for years ite 
30. The thermal properties used in the thermal numerical 
analysis were taken from Boga et al. (1980) and are given in 
Papier 65. 1% 


The transient temperature variation was determined 


based on the following assumptions: 


(i) the thermal properties given in Table 6.1 did not vary 


Significantly with time 
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Table 6.1 Thermal Properties Used in Numerical Analysis 


Thermal Density Specific Heat 
Conductivity 
(k) (p) es) 
Ko /nie nm) ok kg/m? kJ/kg °K 
(Btusznmertoe’R) (ips / Et? )) (Bem 1 be? Ry) 
Crushed Te 209 21625..4 0.638 
Rock/Sand Gt. 17) (S350) (0.20) 
Overburden 6. 230 2082.3 ees) 
efiil) (1.00) (13.050) (OQ, 28) 
oii Sand 64230 2162.4 Tito! 
G7. 00) ins 5 0)) (0.277) 
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Figure 6.1 Finite Element Mesh Utilized for Thermal Analyses 
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ri} the thermal properties were independent of the 
direction of heat flow 

(iii) the ambient air temperature was constant at 10 °C 

(iv) the ambient oil sand temperature was constant at 5 °C 
wy) fluideoflow an the strata,did not’ influence the flow of 


heat 


The results of the thermal analysis are shown in Figures 6.2 
through 6.11. The figures are representative of the 
temperature distribution into the depth of the strata 
uncerlying: the heated foundation at times: 0.25, 7 0i5,% 0.75, 


i, mo, 15, 20, 25 and 30 (years as predicted by (-bni2-7.: 


6.3 Drainage 

The state of drainage in the oil sand strata Baklea the . 
heated foundation was calculated using computer programs: 
VARTEST5 (Appendix E), developed by the author for this 
Study: andieAVINAT, developed by K. J. Bathe “ag Mitt. am 
1 OF Ze 

VARTEST5 calculates the increase in pore fluid pressure 
due to temperature change from the formulations developed in 
Section 3iie3 for each node in the finite element mesh 
employed in the thermal analysis. This program also 
calculates the coefficient of hydraulic GCOUGUCTIVIEY, 
thermal generation rate and the specific storage coefficient 
for each element required as input to the fluid flow program 


ADINAT. The material properties required in the above 
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formulations include: 


e coefficients of thermal volume change for each of the 
components of the oil sand 

e coefficients of compressibility of the components of the 
oil sand 

e dynamic viscosity of the bitumen and water 

®e unit weight (density) of the bitumen, water and sand 


grains 


These values have been determined experimentally or taken 
from the literature and are discussed in detail in Appendix 
A. The step by step operation of VARTEST5 is summarized 
graphically in Figure 6.12. The method in which VARTEST5 
couples the heat transfer program (FEHT2-7) and the fluid 
flow program (ADINAT) is also shown in the figure. 

The fluid flow analysis was performed with a finite 
element mesh consisting of 209 four-node elements. The mesh 
(total of 240 nodes) differs from the heat transfer analysis 
mesh.in that the #edements,.comprising the, crushed. rock 
foundation pad are neglected (the foundation pad is assumed 
to be free draining)2 In® the “analysis the -elements) jwere 
assumed to be linear, that is, having constant permeability 
properties in the iteration process. Two-point Gauss 
numerical integration was employed for the elements. The 
Euler backward method was the numerical time integration 


scheme employed in the transient fluid flow analysis. Bathe 
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Figure 6.12 Computer Program Vartest5 
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(1977) derived the governing equations used with the Euler 
method. 

Two natural boundary conditions were assumed. In the 
first, fluid flow was not permitted normal to the surface 
represented by the axis of symmetry (center of tank). In the 
second, the ground surface was taken as a surface of free 
seepage. Initial conditions were assumed hydrostatic with 
the water table coinciding with the ground surface. 

Under actual field conditions the oil sand mass beneath 
the heated foundation is experiencing simultaneous heating 
and drainage. To ideally simulate these conditions the heat 
transfer and fluid flow program Sou or be excecuted 
frequently using a very small time increment. To accomplish 
this, however, requires an extremely large computational 
effort at unrealistic expense. To provide an optimal balance 
between computational cost and analysis ACCUrACY) @Ehe sort 
Sand strata below the heated foundation wasS permitted to 
heat and drain in one month cycles up to the first year and 
in one year cycles from one year to 30 years. 

A heat and drain cycle begins with the coupling program 
(VARTEST5) reading the elevation heads for each node in the 
finite element mesh. The program then receives the new 
temperature distribution in the Oe sand foundation 
reflecting e"the® propagation. of heat trom ™thewgover tying 
Beetcturesfor™the wtime increment used tor "that" particular 
cycle. The program subsequently reads the previous 


temperature and total head for each node representing the 
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distributions at the end of the prior heat and drain cycle. 
For the first cycle the previous temperature of all nodes is 
> °C-and the pressure distribution is hydrostatic. The new 
temperature distribution is then stored for use in the next 
heat and drain cycle. 

The coupling program then makes a series of 
calculations which include calculating for each node: the 
change in temperature; incremental coefficients of thermal 
expansion for water, bitumen, soil structure and _ sand 
grains; coefficient of compressibility of water, bitumen and 
soil structure; cummulative volume change of water and 
bitumen at new and previous temperature; unit weight of 
water and bitumen at new and previous temperature; change in 
unit weight due to temperature change for water and bitumen; 
in situ viscosity of bitumen, water and fluid; in situ 
hydraulic conductivity (permeability); change in pore fluid 
pressure; change in pressure head; new total head; thermal 
generation rate; and specific storage coefficient. The 
program then calculates the hydraulic conductivity, thermal 
generation rate and specific storage coefficient for each 
element in the mesh. 

The total head, hydraulic conductivity, thermal 
generation rate and specific storage coefficient 
distributions are fed into the fluid flow program which 
determines the total head distribution after the time 
increment of the cycle. The resulting total head 


distribution is converted to excess pore pressure. 
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The results of the analysis are shown in Figures 6.13 
through 6.22, expressed in terms of excess pore pressure 
(pressure exceeding original hydrostatic conditions) with 
Peccnewac cames: 0225, 0.5,10.75, 1,15, tOme to eacO,me> andmao 
years. Figures 6.13 through 6.22 show a band of Sense pore 
pressure propagating into the oil sand strata beneath the 
heated foundation. The shape of these pore pressure bands 
roughly mimick the shape of the isotherms shown in Figures 
6.2 through 6.11. The bands of excess pore pressure in 
figures 6.13 through 6.22 appear to be symetrical beneath 
the tank (diameter of 40 meters) with a region of high pore 
pressure in the center bounded by closely spaced lines of 
equal pore pressure. The center of the excess pore pressure 
bands represent zones of oil sand where drainage is slow and 
the temperature increase is causing pore fluid pressure 
increases. Above the center of the bands the viscosity of 
the bitumen is sufficiently low enough to make it mobile and 
to permit drainage to occur. 

The presence of two thin regions in which the pressure 
Gradients are extremely high may be accounted for in the 
following manner. The high pressure gradient above the bands 
Pewedues in. part to the initial rapid changeso@ hydrauic 
conductivity of oil sand with temperature. Referring to 
Figure A.8 in Appendix A, it is observed that theshydraulic 
Pomduictivitywofeco. la sand changes ,~by. over 3) Gorders of 
magnitude over the temperature range of 5 °C to 25 °C. Below 


the center of the bands of high pore pressure the high 
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pressure gradients reflect the sensitivity of the thermal 
pore pressure generation term with temperature. 

The excess pore presSure distribution with time show 
Significant drainage of the oil sand mass is occurring. Even 
during the time interval of one to 30 years when the length 
of the heat and drain cycles are large (resulting in 
Significant amounts of pore pressure being generated due to 
heating) a large porportion of the oil sand strata beneath 
the tank is almost fully drained. After 5 years the oil sand 
beneath the heated foundation is fully drained to a depth of 
over 20 meters and after 30 years to a depth exceeding 50 
meters (tank center). If the time of the heat aAdedr ain 
cycles were reduced to the one month increment used for’ the 
first) syear, “the porportion of the oil sand mass which is 
drained would be expected to be even greater. The results of 
the analysis indicate that although the rate of pore fluid 
fllowers shower thar tthatitof Theat! flow flit tstesufficient ote 
provide a significant amount of drainage. 

Figures: 603 throtigh #6 (22:alsoeshowsthatethes bands ros 
excess pore pressure in the oil sand strata are generally 
horizontal with respect to the base of the tank (diameter of 
40 meters). This would appear to indicate the heaving due to 
excess pore fluid pressure is fairly uniform and severe 


differential heave is not occurring. 
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6.4 Volume Change 

Volume changes occurring in the oil sand strata was 
assumed to produce only vertical movement (heave) beneath 
the heated structure (Chapter 4). The magnitude of heave was 
calculated by considering volume changes resulting from 
volume change of pore fluid in response to relief of pore 
fluid pressure which exceeds the total vertical confining 
stress (weight of structure and overlying sediment) and 
thermal volume change of the sand grains and soil structure. 

The pore fluid pressure distribution was calculated 
twice for each increment of time, after the heating cycle 
and after the draining cycle. During the heating cycle the 
pore fluid pressure response was determined assuming 
undrained conditions with no volume change using equation 
sozayer' wThe tporesefluad dpressure) distrabution¢duridng the 
drainage cycle was predicted by the fluid flow analysis. In 
either case, if the pore fluid pressure exceeded the total 
vertical confining pressure the strata was permitted to 
expand such that the pore fluid pressure balanced the total 
stress. This volumetric expansion was converted to heave by 
assuming no lateral strain. If during subsequent time 


increments, zones of pore fluid pressure exceeding the total 


‘An inconsistancy 1s apparent Since calculation of pore 
fluid pressure from equation 3.24 ignores the volume change 
due to the thermal expansion of sand grains and soil 
structure. This error, however, is slight since the thermal 
volume change of sand grains and soil structure is 
relatively small (less than 1% at 175 °C, Figure 5.11) and 
thus the bitumen and water porosities (ratio of volume of 
bitumen or water to volume of oil sand mass) in equation 
3.24 would change only slightly. 
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vertical confining stress became sufficiently permeable to 
permit full drainage of pore fluid, the heave associated 
with this zone of pore pressure was neglected (i.e., the 
soil structure was allowed to collapse). In other words, the 
heave associated with pore fluid pressure relief was 
determined incrementally with each heating and draining 
cycle. A plot of heave beneath the heated tank resulting 
from pore fluid pressure relief is shown in Figure 6.23. 
Figure 6.23 indicates that heave resulting from pore 
fluid pressure relief is sensitive to the time increment 
employed in the analysis. For the first year when the heat 
and drain cycles were limited to one month, a maximum naaes 
in the order of 0.7 cm (tank center) was observed at 3 
months and steadily decreased with time to one year. 
However, when the heat and drain cycles were increased to 
one year a heave of about 3 cm (tank center) was found after 
5 years and decreased with time to approximately 0.1 cm 
after 30 years. This dependence on the time increment is due 
almost entirely to the amount of pore fluid pressure 
Calculated during the heating portion of the heat’ and drain 
Cycle. Examination of equation 3.24 indicates that the pore 
fluid pressure ‘change, Au, is dependent ‘on the) temperature 
ehange,) AT. During the |one year, “Gime increment the 
temperature change is much greater than for a one month time 
increment and, as a result, significantly more pore fluid 
pressure is generated. One | would’ | expect (apt sche  itime 


increment was held constant throughout the analysis at one 
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month, the amount of heave resulting from pore fluid 
pressure relief would decline steadily throughout the 
analysis. The maximum differential heave resulting from pore 
fluid pressure relief was found to be about 1.4 cm and 
occurred after 5 years of heating. 

The heave associated with thermal volume change of the 
Sand grains and soil structure with time is shown in Figure 
6.24. A maximum heave of about 35 cm (tank center) was 
observed after 30 years of heating. Unlike the heave 
associated with pore fluid pressure relief, the heave 
increases at any point under the tank progressively with 
time and appears far less dependent on the time increment. 
This is because the thermal volume change of sand grains and 
soil structure depends only on the volume of oil sand 
experiencing a temperature increase and 1s independent of 
the amount of pore fluid drainage occurring. It 1S apparent 
that the heave associated with thermal volume change of the 
Sand grains and soil structure is one order of magnitude 
Greater’ than ithat due to pore fluid presstre reliet. The 
Maximum differential heave resulting from thermal volume 
change of the sand grains and soil structure was found to be 
about 14 cm and occurred after 30 years of heating. 

Figure 6.25 shows the total heave beneath the heated 
foundation with time. Comparison of this plot with Figures 
6.23 and 6.24 shows the vast majority of heave is due to the 
thermal volume change of the sand grains and soil structure. 


A maximum differential heave of 14 cm was observed after 30 
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yearsimof Bheadting!, Thisseindicatesrepores filmid pressure 
gevelopede during! lheatingonofl then oidmeisandsstrataadrains 
sufficiently quickly enough to negate large volume changes 


associated with the undrained condition. 


6.5 Conclusions 

The analytical model presented in this study 
investigated the effect of heat flow from an 80 m diameter 
Storage tank on underlying oil sand strata. The tank was 
assumed to rest on a conventional crushed rock and_ sand 
foundation pad, 2 m in thickness, and held at a constant 
bemperature of) ot? 50 n*# Cees Thesuresuitsrcofvothes. theoretical 
analysis showed a band of excess pore pressure (pressure 
exceeding original hydrostatic conditions) propagating into 
the oil sand strata beneath the heated foundation. The 
center of the excess pore pressure bands represent zones of 
oil sand where drainage is slow and the temperature increase 
is causing pore fluid pressure increases. Above the center 
of the bands the viscosity of the bitumen is sufficiently 
low enough to make it mobile and to permit drainage to 
eceun.. 

The excess pore pressure distribution with time show 
Significant drainage of the» oil» sand mass 15 occurring, 
After 5eyears ofsheatingatheroil ysand beneathertheedheated 
foundabionovisimitily drained! torasdepthnrol over Zopmetersy 


andtafteres0iyearseto am depthimexceedingsi50tometers atitank 
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The bands of excess pore pressure in the oil sand 
Sstratanare)generally horizontal .with»respect),to+the ~base, of 
the tank indicating the heaving due to excess pore fluid 
pressure is fairly uniform and severe differential heave is 
noOEfoccurring, 

The results of the analysis indicate that although the 
rate of pore fluid flow is slower. than that of heat Plow Lt 
is sufficient to provide a significant amount of drainage. 

Volume changes occurring in the oil sand strata was 
assumed to produce vertical movement (heave) beneath the 
heated structure. The magnitude of heave was calculated by 
considering volume changes resulting from volume change of 
pore fluid in response to relief of pore fluid pressure 
which exceeds the total vertical confining stress (weight of 
structure and overlying sediment) and thermal volume change 
of the sand grains and soil structure. 

The heave resulting from pore fluid pressure relief was 
found sensitive to the time increment employed in the 
analysis. For the first year when the heat and drain cycles 
were limited to one month, a heave in the order of 0.7 cm 
(tank center) was observed at 3 months and steadily 
decreased with time to one year. When the heat and drain 
cycles were increased to one year, a heave of about 3 cm 
(tank center) was found after 5 years and decreased with 
time to approximately 0.1 cm after 30 years. This dependence 
on the time increment is due almost entirely to the amount 


of pore fluid pressure generated during the heating portion 
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of the heat and drain cycle. If the time increment was held 
constant throughout the analysis at one month, the amount of 
héeaveseresulting! from® poressfbuidiupréssure® feéLivet Swotlld 
decline® steadily throughout “the! analysis? The maximum 
differential heave resulting from pore fluid pressure relief 
was found to be about 1.4 cm and occurred after 5 years of 
heating. 

The heave associated with thermal volume change of the 
sand grains and soil structure with time showed a maximum 
heave of about 35 cm (tank center) occurred after 30 years 
of heating. Unlike .the heave associated with pore fluid 
pressure détiveg tute heave increases at any point under the 
tank progressively with time and is far less dependent on 
the time increment. The heave associated with thermal volume 
change of sand grains and soil structure is one order of 
Magnitude greater than that due to pore fluid pressure 
relief. The maximum differential heave resulting from 
thermal volume change of the sand grains and soil structure 
was found to be approximately 14 cm and occurred after 30 
years of heating. 

When considering the total heave beneath the heated 
foundation with time, the vast majority of heave is due to 
the thermal volume change of the sand grains and soil 
Structure. A maximum differential heave of 14 cm was found 
after £30 diyearal to fe heating Ne This tindicaves Stpore Veihund 
pressure generated during heating of the oil sand strata 


drains fast enough to negate large volume changes associated 
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with the undrained condition. 

The theoretical analyses used in this study did not 
include the influence of gas exsolution. The amount of gas 
driven out ‘of solution from isthe fluid phase ofsoil ssandy is 
dependent on the kind and amount of gas in solution, the 
increase in temperature and the pore fluid pressure. The 
amount of gas initially in solution in the fluid phase of 
the oil sand is dependent on the initial pore pressure. 
Since hydrostatic pore pressure conditions initially exist 
in the surface mineable area little gas is expected to be in 
solution near the ground surface but the amount increases 
progressively wit depth. The increase in total stress 
exerted by the weight of the structure and foundation pad 
reduces the pore volume increase from an increase in 
temperature and therefore has a significant effect on the 
amount of gas driven from solution. Since the propagation of 
heat from the structure is extremely slow when compared to 
that experienced in the laboratory, it would be difficult to 
experimentally assess the amount of gas drainage which would 
occur under actual field conditions. The permeability of oil 
sand to gas would also be difficult to assess but may be 
Significant. It would appear that some amount of volume 
change due to gas exsolution could be expected to contribute 
to heaving beneath the heated structure, the amount of which 
would be difficult to determine. Since the isotherms beneath 
the tank were shown to_ be somewhat horizontal the 


contribution to differential heave may not prove to be of 
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major concern. For the above reasons, it was decided to not 
includempgas exsolution)! inputhisi=s first’ analysis off)the 
performance of heated foundations on oil sand. 

Variations in Stratigraphy beneath the heated structure 
may have a significant influence on the amount of heave 
which occurs. A simplified stratigraphy consisting of a thin 
Pill@ layer) overlying’ rich oil sand was’assumed for this 
Study. This stratigraphy was selected because rich oil sand 
which contains a large percentage of bitumen has an 
extremely low permeability under formation temperatures and 
therefore when subjected to an increase in temperature would 
generate large amounts of pore pressures and possibly 
represent a worse case. It was of concern to assess the rate 
at which these high pore pressures would dissipate. The 
results of the analyses shows a significant portion of this 
pore fluid pressure was dissipating and the heave associated 
with pore fluid volume increase was not significant. 

Should clay shale layers be present in the oil sand, 
drainage would initwally be impeded, however, the 
permeability of the clay shale layer would possibly increase 
dramatically through fracturing: due top ernermals volume 
Changes. off the underlying  o11 Sand.” Inejcontrast athe 
presence of lean or bitumen free sand would Significantly 
improve the hydraulic conductivity of the strata since the 
dynamic viscosity of water is: much’ Tower? sthan that of 
bitumen. In addition, the- amount of pore fluid pressure 


generated in lean oil sand would be less than for rich oil 
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sand. 

Complex stratigraphy such as discontinuous layers 
adversely affect the amount of differental heave which would 
occur and would have to be modelled for the specific 


conditions. 
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7. SUMMARY AND CONCLUSIONS 


7.1 Summary 

me yextractionsand jsupgrading sof, jybitumenw ata sumtace 
mining projects may require the construction of heat 
generating structures above oil sand strata. To date, the 
design procedure adopted involves limiting the temperature 
in the oil sand beneath the structure to 45 °C. The design 
criteria was based on undrained triaxial testing of 
disturbed oil sand over the temperature range of room 
temperature, toe 904 °C. ,For astructures. generating: .large 
amounts of heat, such as heated bitumen storage tanks, the 
design criteria results in very complex and expensive 
foundation pads utilizing ventilation systems. 

The purpose of this study was to provide a better 
understanding of the geotechnical behavior of oil sands 
under conditions that develop after the placement of heated 
structures on oil sand and to provide a rational basis of 
design for foundation pads for these structures. 

The analytical concepts presented in this study provide 
the rationale to assess the influence of heat on oil sand 
(thermal volume change and pore fluid pressure response), 
rate Tof -heat+eflowneand -ratecandyextentyofaftluid: Elow: (or 
drainage). These concepts give the basis for the theoretical 
analysis of the geotechnical behavior of oil sand underlying 
heated foundations. The- theoretical analysis serves to 


couple heat flow~theory to the fluid flow expression to 
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determine whether the oil enna masS will be fully drained, 
fully undrained, or partially drained. If partial drainage 
occurs, the pore fluid pressure change and volume change can 
be calculated. The theoretical analysis can be further 
refined by considering the influence of gas exsolution and 
pore fluid vaporiztion or distillation and by considering 
the three dimensional case. 

The material properties required for the analyses are 
coefficients of thermal volume change and the coefficients 
of compressibility for each of the components of oil sand, 
heat flow parameters, as well as the dynamic viscosity of 
the bitumen and water and the density of bitumen, water and 
sand grains. These values have been determined 
experimentally or taken from the literature. 

The analytical model presented in this study 
investigated the effect of heat flow from an 80 m diameter 
storage tank on underlying oil sand strata. The tank was 
assumed to rest on a conventional crushed rock and sand 
founda tonirpad / {2 imei m:vhroknessiftand) Cheld t*at Sta "constant 
temperature: .of i750 °Ge THe tgeneral© stratigraphy "chosen 
Gonscisitied lo fiieheecqranilar Sipad Sromerrop., of 3 miof. tilt 
overlying oil sand. 

The results found in this study are based on a specific 
Stratigraphy and measured properties for one oil sand. If 
different stratigraphy and properties of oil sand, such as 
significant gas exsolution, exist then those particular 


conditions would- have to be modelled. 
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7.2 Conclusions 

The results of this study showed that when examining 
the foundation performance beneath a heated structure, heave 
and differential heave is the sole factor affecting the 
foundation performance of the oil sand. The possibility of 
shear failure and lateral movement of the oil sand and 
structure are extremely remote. Settlement of the structure 
does not take place even under successive heating and 
cooling cycles. 

An expression was developed to calculate the net volume 
change of an oil sand mass subjected to a change in 
temperature under fully drained conditions. The volume 
change of the oil sand mass arises from the thermal 
expansion of the individual sand grains and structural 
changes in the sand skeleton. 
| An equation which allows the volume change of oil sand 
Subjected to increases in temperatures under completely 
Taeeaned conditions to be calculated was presented. The 
volume change of the oil sand mass is due to expansion of 
the individual sand grains and pore fluids as well as 
structural changes in the sand skeleton. 

An expression was derived which expresses the volume of 
pore fluid which would drain from a fully saturated oil sand 
mass experiencing a change in temperature. The amount of 
fluid expelled from an oil sand mass during a drained 
thermal expansion test -is the sum of the volume change of 


the pore fluids less the change in void volume and 
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structural volume changes. 

An expression was developed to calculate the increase 
in pore fluid pressure in a fully saturated oil sand mass 
subjected to a temperature change under fully undrained 
conditions. 

The theoretical expressions to evaluate the rate of 
fluid from heating were derived and are presented. The 
method utilized to solve the governing differential fluid 
flow equation for two dimensional transient fluid flow 
beneath a heated foundation using the finite element 
technique is described. 

The theoretical formulations to determine the rate of 
heat flow were derived and are presented. The solution to 
the heat equation for two dimensional transient heat flow 
beneath a heated foundation utilizing the finite difference 
technique is illustrated. | 

The expressions developed in this study provide a 
method to analyze the rate at which heat flows through an 
oil sands mass, the rate at which excess pore fluid pressure 
resulting from an increase in temperature dissipates, and 
the amount of volume change which occurs. 

Equipment was successfully designed and constructed 
with the capability to explore the one-dimensional thermal 
expansion’ «of oil sands \andeits constituent components to 
Cempetatifes approaching 200m@°Ceandtpressuresoine “excescsn Wor 


3.5 MPa. A brief summary of each constituent of the 


equipment is presented. 
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The test procedures which were developed for use in the 
laboratory program worked successfully and are described. 
The method utilized to prepare oil sand samples from frozen 
core resulted in very high quality test specimens. The 
procedure adopted is illustrated. 

Thermal expansion tests performed on water resulted in 
volume increases of about 15 per cent at 200 °C. The thermal 
volume change of water did not vary significantly for the 
pressure levels employed. Excellent agreement between 
experimental and published data indicates the pressure cell 
employed was calibrated and operating properly. 

Thermal expansion experiments carried out on hot water 
extracted bitumen resulted in volume increases of about 15 
per cent at 200 °C. The test results indicate that the 
Be eeiavive volume change of bitumen over the temperature 
range considered is not pressure dependent. The volume 
change of this dead bitumen is more linear with temperature 
than that of water. 

For each thermal expansion test conducted on compacted 
tailing sand, a negative volume change was observed for the 
first heating increment (room temperature to 50 °C). The 
loosest sand sample (initial density of 1.46 Mg/m* at room 
temperature) gave the greatest volume decrease (0.36%) over 
this temperature range. The sample initially compacted to a 
density of 1.62 Mg/m*? at room temperature showed a volume 
decrease of 0.08% over- the temperature range from room 


temperature to 50 °C. After heating to 200 °C, this sample 
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was cooled to room temperature and a volume decrease of 
0.10% was measured indicating that the volume decrease was 
due to a compaction of the soil structure. Relative movement 
Or reorientation of sand grains must have occurred during 
the temperature increase to 50 °C to allow an overall volume 
decrease. Over the temperature range of 50 to 200 °C all 
samples tested displayed a volume increase roughly paralling 
what ®ofP%%a “quartz "erystal, b'indicating*thecsoil structure 
achieved a tight packing arrangement at 50 °C and the 
Subsequent volume increase with temperature results only 
from the thermal expansion of the individual grains. 

Undrained heating of oil sand to 200 °C will result in 
a volume increase of about 6 percent. Most of this increase, 
oer 5 Ypercenit, S-1swedue = foi. the Spore fluids, water and 
bitumen, with the bitumen expanding slightly more than 
water. The remainder of the volume increase is due to the 
thermal expansion of the sand grains. The dense, 
impenetrative structure of in situ oil sand will show a 
negligible amount of thermal volume change by itself. 

From the tests performed on oil sand the thermal volume 
increase appears to be independent of the effective 
confining pressure for the pressure range considered in this 
Study. For the oil sand tested here, significant thermal 
volume increases would be expected under back pressures 


below 2000 kPa because of gasS exsolution at temperatures 


above 180 °C. 
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Results obtained from the undrained and drained thermal 
expansion tests performed on oil sand were used to generate 
a plot representative of the thermal expansion of in situ 
bitumen. The results indicate that in situ gas _ saturated 
bitumen may show temperature volume increases 50% greater 
than the extracted gas free bitumen. 

The results of the theoretical analysis showed a band 
of excess pore pressure (pressure exceeding original 
hydrostatic conditions) propagating into the oil sand strata 
beneath the heated foundation. The center of the excess pore 
pressure bands represent zones of oil sand where drainage is 
Slow and the temperature increase is lceneane pore ffluad 
pressure increases. Above the center of the bands the 
viscosity of the bitumen is sufficiently low enough to make 
it mobile and to permit drainage to occur. 

The excess pore pressure distribution with time show 
Significant drainage of the oil sand mass is occurring. 
After 5 years of heating the o11 sand beneath the heated 
foundation is *fully ‘drained tova depth of over’ 20 meters, 
and after 30 years to a depth exceeding 50 meters (tank 
center). 

The bands of excess pore pressure in the oil sand 
strata are generally horizontal with respect to the base of 
thestank indicating the heavingsdte Oto Mexcesseeporerviluta 
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The results of the analysis indicate that although the 
paterofipore fluid flow is slower than that of heat ‘flow, it 
is sufficient to provide a significant amount of drainage. 

Volumes changeso\ocecurringrein(| tbheionl ssanddstrata ivas 
assumed to produce vertical movement (heave) beneath the 
heated structure. The magnitude of heave was calculated by 
considering volume changes resulting from volume change of 
pore fluid in response to relief of pore fluid pressure 
which exceeds the total vertical confining stress (weight of 
Structure and overlying sediment) and thermal volume change 
of the sand grains and soil structure. 

The heave resulting from pore fluid pressure relief was 
found sensitive to the time increment employed in the 
analysis. For the first year when the heat and drain cycles 
were limited to one month, a heave in the order of 0.7 cm 
(tank center) was observed at 3 months and _ steadily 
decreased with time to one year. When the heat and drain 
cycles were increased to one year, a heave of about 3 cm 
(tank center) was found after 5 years and decreased with 
time to approximately 0.1 cm after 30 years. This dependence 
on the time increment is due almost entirely to the amount 
of pore fluid pressure generated during the heating portion 
of the heat and drain cycle. If the time increment was held 
constant throughout the analysis at one month, the amount of 
heave resulting from _ pore fluid pressure relief would 
decline steadily throughout the analysis. The maximum 


differential heave resulting from pore fluid pressure relief 
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was found to be about 1.4 cm and occurred after 5 years of 
heating. 

The heave associated with thermal volume change of the 
sand grains and soil structure with time showed a maximum 
heave of about 35 cm (tank center) occurred after 30 years 
of heating. Unlike the heave associated with pore fluid 
pressure relief, the heave increases at any point under the 
tank progressively with time and is far less dependent on 
the time increment. The heave associated with thermal volume 
change of sand grains and soil structure is one order of 
magnitude greater than that due to pore fluid pressure 
relief. The maximum differential heave resulting from 
thermal volume change of the sand grains and soil structure 
was found to be approximately 14 cm and occurred after 30 
years of heating. 

When considering the total heave beneath the heated 
foundation with time, the vast majority of heave is due to 
the thermal volume change of the sand grains and soil 
structure. A maximum differential heave of 14 cm was’ found 
after 30 years of heating. This indicates pore fluid 
pressure generated during heating of the oil sand _ strata 
drains fast enough to negate large volume changes associated 
with the undrained condition. 

The results of the study showed significant drainage of 
the oil sand mass was occurring and that the vast majority 
of vertical surface movement or heave was due to the thermal 


expansion of the sand grains and soil structure. The amount 
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of differential heave predicted from the theoretical 
analysis appears manageable when incorporated into a 


Structural design of a storage tank. 


7.3 Recommendations for Future Work 

The oil sand used in the experimental testing was 
obtained from core samples taken from an oil sand outcrop at 
Saline Creek, approximately 1 km south of Fort McMurray. It 
is apparent that the stresses in this outcrop of oil sand 
have been reduced at a very slow rate during geologic time 
and therefore, significant portion of the gas initially 
present in the oil sand has drained. When considering the 
design of a heated structure on oil sand, testing of the oil 
Sand should be conducted on a site specific basis to 
determine if similar amount of gas drainage has occurred. It 
is possible with the apparatus used in this testing program 
to monitor and sample gases from oil sand tested at high 
temperatures. 

Experimental testing should be undertaken to provide a 
better understanding of the heat flow properties of oil sand 
ana to. measure the change «in hydraulicucondtctivity of o1) 
Sand with temperature. 

Should strata other than oil sand be present in the 
foundation beneath the heated structure, similar testing on 
these soils should be conducted. 

The theoretical analysis used in this study can be 


further refined by considering the influence of gas 
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exsolution and pore fluid vapourization or distillation and 
by considering the three dimensional case. The influence of 


closely Spaced heated structures should also be 


investigated. 
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APPENDIX A - THEORETICAL DEVELOPMENT 


Introduction 

The engineering properties of the oil sand constituent 
components which were required as input for the numerical 
analyses were determined either from the experimental 
testing OF from published information. This appendix 
outlines briefly the source of the data used to derive each 
property and how each property varies with temperature and 


pressure. 


Incremental Coefficient of Thermal Expansion 


(i) Water 

The incremental coefficient of thermal expansion for 
water was derived from both experimental and published data. 
The coefficient, calculated based on 10 °C temperature 
increments, did not very Significantly over the pressure 
range of 0 to 3500 kPa (temperature range of 5 to 200 °C). 
As a result, it was possible to derive a unique relationship 
between the incremental coefficient of thermal expansion and 
temperature. This relationship is shown graphically in 


Figure A.1. 
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The incremental coefficient of thermal expansion for in 
Situ bitumen was derived solely on data collected during the 
experimental testing. The coefficient, calculated based on 
10 °C temperature increments, is assumed to be pressure 
independent over the temperature range considered in this 
study (5 to 200 °C). A plot of the: incremental coefficient 
of thermal expansion with increasing temperature is shown in 


Figure A.1. 


Melty) tin SitweSothsStructure 

The relationship between the incremental coefficient of 
thermal expansion for the in situ oil sand soil structure 
and temperature iS given in Figure A.1. The coefficient was 
calculated considering data obtained from drained thermal 
expansion tests performed on high quality undisturbed oil 
Sand samples and adapting these results to approximate in 
Situ conditions through density extrapolation (Chapter 5). 
The incremental coefficient of thermal expansion for the 


Soll structure was calculated based on temperature 


menementsaof. i2n5mec. 


(iv) Sand Grains 


The incremental coefficient of thermal expansion for 
the oil sand grains was calculated based on data extracted 
from the 1966 edition of the Handbook of Physical Constants 
published by the Geological Society of America, Inc. A plot 


of the incremental coefficient of thermal expansion with 
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increasing temperature is shown in Figure A.1. 
Incremental Coefficient of Compressibility 


(i) Water 

The coefficient of compressibility for water was 
calculated based on data taken from the 1968 edition of the 
Japanese Society of Mechanical Engineers Steam Tables. Over 
the pressure range of 0 to 3500 kPa the coefficient did not 
weey for temperature levels-up to and) including: ~200 °G@. 
Therefore, it was possible to derive a unique relationship 
between the incremental coefficient of Panprees ee and 
temperature. This relationship 1S shown graphically on 


Figure A.2. 


(}i) In Situ Bitumen 

The incremental coefficient of compressibility for in 
Situ bitumen was calculated based on information extracted 
from the work done by Svrcek and Mehrotra (1981). In situ 
bitumen was assumed to be gas saturated with 80 percent 
methane and 20 percent carbon dioxide. As with the case of 
water, over the pressure range under consideration the 
coefficient is not pressure sensitive and as a result it was 
possible to derive a unique relationship between the 
incremental coefficient of compressibility and temperature 
(Figure A.2). Since only. data to 100 °C was available, the 


curve was extrapolated to 200 °C as shown in the figure. 
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(iil) Soil Structure 

The incremental coefficient of compressibility for the 
in situ 011 sand soil structure was assumed to be both 
pressure and temperature independent for the pressure and 
temperature ranges utilized for this study. The value used 
in the analysis was selected from the range suggested by 


Dusseault (1980) for in situ oil sand (Figure A.2). 
Dynamic Viscosity 


(i) In Situ Bitumen 

In situ bitumen, in this study, 1s assumed to obtain 80 
percent methane and 20 percent carbon dioxide in solution. 
Dynamic viscosity-temperature relationships for methane and 
carbon dioxide saturated bitumens, for various pressure 
levels, have been reported by Svrcek (1979). 

Plots of the dynamic viscosity of methane saturated 
bitumen over the pressure range of atmospheric (gas free 
oil) to 2.5 MPa is shown in Figure A.3. Over this pressure 
range the dynamic viscosity of methane saturated bitumen is 
relatively insensitive tor pressure changes. As a result, it 
waS posSible to derive a unique relationship between the 
dynamic viscosity of methane saturated bitumen and 
temperature (Figure A.3) without introducing any appreciable 
amount of error. 

Plots of the dynamic viscosity of carbon dioxide 


Saturated bitumen over the pressure range of atmospheric 
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(gas free oil) to 2.34 MPa is shown in Figure A.4. Over this 
pressure range the dynamic viscosity of carbon dioxide 
Saturated bitumen is significantly pressure dependent over 
ene temperature interval of 5 to 75 °C: Since: the! (majority 
of pressure dependency occurs in the room temperature range 
(where the bitumen id known to be virtually immobile) and 
considering that in situ bitumen is assumed to contain only 
20 percent carbon dioxide, it iS reasonable to assume a 
unique relationship between the dynamic viscosity of carbon 
dioxide saturated bitumen and temperature as shown in Figure 
A.4. 

Figure A.5 gives the dynamic viscosity-temperature 
relationship for in situ bitumen based on the assumption 
that in situ bitumen contains 80 percent methane and 20 


percent carbon dioxide in solution. 


(if) Water 

The dynamic viscosity of water with increasing 
temperature was obtained from the 1980 edition of 
Thermodynamics and Transport Properties of Fluids by Rogers 
and Mayhew. A plot of -the. dynamic viscosity-temperature 


relationship for water is*given in Figure A.6. 


(iii) Fluid (Bitumen and Water) 
The dynamic viscosity-temperature relationship for the 
fluid was based on an oil -sand having 15 weight percent 


bitumen (31 volume percent) and 2 weight percent water (4 
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volume percent). This relationship is shown graphically in 


Figure A.6. 
Unit Weight (Density) 


(i) In Situ Bitumen 

The room temperature unit weight for bitumen containing 
80 percent methane and 20 percent carbon dioxide was assumed 
to be 10.104 kN/m? (1.03 g/cm?) (Svrcek and Mehrotra, 1981). 
Utilizing the volume change-temperature relationship derived 
in the laboratory testing (Chapter 5) a curve representing 
the change in unit weight with temperature for in situ 
bitumen was calculated. This relationship is presented in 
Figure A.7. 

Based on the data reported by Svrcek ane Mehrotra 
(1981) the presence of carbon dioxide in solution has little 
effect on density over the pressure range considered in this 
Study. Similarly, the unit weight of methane’ saturated 
bitumen does not vary significantly over the pressure range 
used in this study (Svrcek and Mehrotra, 1981). Therefore, 
the unit weight-temperature relationship given in Figure A.7 


igs valid for analytical purposes. 


(ji) Water 


The unit weight-temperature relationship for water was 
derived from both experimental and published data. Since the 


volume change of water did not vary. significantly over the 
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pressure range of 0 to 3500 kPa (temperature range of 5 to 
200 °C) it was possible to derive a unique relationship 
between unit weight and temperature. This relationship is 


shown graphically in Figure A.7. 


(fii) Fluid (Bitumen and Water) 

The unit weight-temperature relationship for the fluid 
was based on an oil sand having a 15 weight percent bitumen 
(31 volume percent) and 2 weight percent water (4 volume 
percent). This relationship is shown graphically in Figure 


Becks 


Hydraulic Conductivity 

The hydraulic conductivity- temperature relationship 
for the fluid (bitumen and water) was derived from the unit 
weight-temperature and dynamic viscosity-temperature 
relationships and assuming an absolute permeability of 3 x 
fos '* me (Chapter. 3). A plot “of the ~hydraulic # hydraudie 


conductivity with increasing temperature is shown in Figure 


A.8. 
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APPENDIX B - DETAILED DESCRIPTION OF APPARATUS 


Consolidometer 


(i) Sample Jacket 

The consolidometer sample jacket is composed of 420 
series Stainless steel with a coefficient of linear thermal 
expansion of 10.8 micrometers per meter per degrees celsius 
[Metals Handbook, Vol. 3, 9th Edition). The jacket was 
tempered to increase its hardness in order to. prevent 
scratching of the metal by the sand grains. Also made 
smooth. 

The sample jacket has an inside diameter of 7.63 cm and 
a wall thickness of 0.635 cm. The jacket was designed to 


accommodate specimens ranging in height from 2.0 to 4.0 cm. 


(ii) Consolidometer Piston 

The consolidometer piston consists of 304 series 
Stainless steel with a coefficient of jJianear thermal 
expansion of 17.8 micrometers per meter per degrees celsius 
(Metals Handbook, Vol. 3, 9th Edition). 

The clearance between the piston and sample jacket is 
approximately 0.025 cm. The overall weight of the 
consolidometer piston (including concrete cylinder) is 38.84 


N. Location and dimensions of thermocouple and drainage 


ports are given in Appendix D. 
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There are two closely spaced sealing O rings located on 
the lower portion of the piston and a rider ring positioned 
on the upper portion. The purpose of the latter is to ensure 
vertical movement of the piston and to reduce frictional 


resistance between the piston and sample jacket. 


(iii) O and Rider Rings 

All O rings used with the apparatus are composed of a 
chemical compound commercially referred to as Viton. Viton, 
unliké conventional o-ring material, is resistant to 
degradation over the operational temperature range of the 
apparatus and in addition is compatible with petroleum 
products. 

Teflon, a chemically stable, heat resistant plastic 
polymer, was selected as the material for the rider ring. 
Teflon proved to be an ideal material owing to its low 
coefficient of friction, good commercial availability and 


ease of machining. 


(iv) Pedestal 


The pedestal is composed of the same material as _ the 
piston and contains 2 drainage ports, location and 


dimensions of which are given in Appendix D. 


(v) Porous Stones 


The consolidometer contains three porous stones. Two 


0.32 cm thick porous stones are imbedded within the piston 
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as shown in Figure 2.1 and one 0.64 cm thick porous stone in 


the pedestal. 


(vi) Filter Paper 

The filter paper used in the consolidometer consists of 
very finely wooven glass fibers. Unlike conventional filter 
paper, glass fiber paper is heat resistant. The filter paper 
has the finest mesh (effective retention, 1.2 x 10°° m) of 
those commercially available giving the greatest amount of 


flow resistance. 
Measuring Devices 


(ij) Pressure Transducers 

The transducers used to measure the back and air 
cylinder diaphram pressures were commercially produced by 
the Viatran Corporation. Briefly, the transducers have 
diaphrams which deflect when subjected to fluid under 
pressure. This deflection is measured by foil strain gauges 
mounted on the diaphram and converted to output voltage 
Signals. The design is such that the output signal varies 
binearly swithewithsapplied pressures 

The air cylinder pressure transducer (Model 222) has an 
operational range of 0 through 2.0 MPa. The back pressure 
transducer (Model 108) has an operational range of 0 to 7.0 


Mpa. Detailed specifications are readily available from the 


manufacturer. 
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(ii) Linearly Varying Displacement Transducer 

The LVDT used in the apparatus was manufactured by the 
Hewlett Packard Company. 

The LVDT performs on the following theory: the direct 
current input is transformed into alternating current by an 
oscillator. The alternating current excites the primary 
windings inducing voltage in the secondary windings, the 
amount of which is determined by the location of the axial 
core. Secondary circuits enable the output to be read as DC 
voltage which is linearly proportional to the axial core 
displacement. 

The model used in the apparatus was the HP 24DCDT. 
Detailed specifications are readily available from. the 


manufacturer. 


(iii) Thermocoup]es 

The thermocouples used in the apparatus were 
manufactured by the Barber-Colman Company. The thermocouples 
are Type vu (iron-constantan element) with fiberglass 
insulation. The element is housed in 0.3175 cm_ stainless 
Stee berprotectives i bubesrwhich{ lus leclosed thaumathemendseThe 
operational temperature range of the thermocouple is from 0 
to approximately 400 °C. Detailed specifications are readily 


available from the manufacturer. 
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The band heater used in the apparatus was obtained from 
Thermel Incorporated. The Thermastrip heater band (Model 
CD-25350) contains two elements and is rated at 650 watts. 
The band has an inside diameter of 8.9 em; is 6.35 cm wide 
and operates on 115 volts. The heater also has a 6.35 cm 


terminal box. 


Diaphram Air Cylinder 

The diaphram air cylinder was obtained from the 
Bellofram Corporation. Its theory of operation was discussed 
in Chapter 4. The following is a brief summary of the 


pertinent specifications: 


TYPES 5 
Si ZE% 24 
SERIES: E 
ROD: BP 
BORE: | ioe 77cm 
STROKE: 6.6CRem 
MAXIMUM APPLIED AIR PRESSURE: 1260 kPa 


Further specifications are available from the manufacturer. 


Temperature Controller 
The temperature controller (Model 580) is 
microprocessor-based and manufactured by the Barber-Colman 


Company. The operational use of the instrument is discussed 


in Chapter 4. 
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Theacontrollertihastea tratedlidectragy ofmpe(0y5% lof 
reading + 2 °C) and a temperature range of 0 to 400 °C. The 
controller has a temperature sampling time of 3 per’ second. 


Detailed specifications are available from the manufacturer. 


Insulating Shell 

The material composing the insulating shell was 
manufactured by the Carborundum Company Insulation Division. 
The Fibrefrax Moist-Pack D used to construct the shell is 
made of ceramic fiber and an inorganic binder which is 
easily workable but~ cures to a rigid structure. The 
insulation is composed primarily of Alumina-Silica with an 
operational limit of 1260 °C. The material has good thermal 
reflectance, low heat storage and low thermal conductivity. 

The insulating shell itself is 1.27 cm thick and fully 
encompasses the consolidometer. Detailed specifications are 


available from the manufacturer. 


Data Aquisition Unit 

The data aquisition unit was manufactured by the 
Hewlett Packard Company. The unit consists of the HP-85F 
Desktop Computer and the HP-3497A Data Aquisition Control 
Unit. 

The data logger is programmable to print measurements 
on paper as well as record on magnetic tape. During testing 
up to 3 channels can be-monitored on the CRT display. The 


data recorded by the unit on the magnetic tape is 
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transferable directly to the University computer for 
processing. 

The operational use, although not complex, is  smewhat 
involved and will not be discussed here. Procedural 
instructions and unit specifications are available in the 


Operator’s Handbook accompanying the equipment. 
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APPENDIX C - LABORATORY CALIBRATIONS 


Calibration of Measuring Devices 

To process the data collected by the data acquisition 
unit it was necessary to calibrate each of the measuring 
devices independent of the apparatus. Each of the measuring 
devices had an operational range far exceeding that utilized 
by the apparatus but, to precisely define the calibration 


functions for each instrument, each was taken to capacity. 


(i) Calibration of Pressure Transducers 

Each pressure transducer was calibrated. using the 
Barnet Instruments Ltd. Calibrator. Documentation concerning 
its use is available with the instrument and will not be 
discussed here. Figure C.1 shows the calibration curve for 
the 7.0 MPa pressure transducer (used to monitor the back 
pressure within the apparatus), while Figure C.2 illustrates 
the curve for the 2.0 Mpa transducer. Since the 2.0 Mpa 
transducer cannot directly measure the consolidometer piston 
pressure but only the applied air pressure to the diaphram 
air cylinder, an additional calibration was required. It 
involved calibrating the applied oe pressure to. the 
diaphram air cylinder with the load exerted on the piston. 
This calibration was performed using a standard proving ring 


and 1S showm in Pigure’ C.3. 


188 


sworranatdaQ YROVANOEAS - D xTOWETGA J ECO T i 


A | P 
eedives pnizeseM te coi tend 
sjsh etd yd Begosliog s7sb reas ‘enex07g oT 
nS or 


a ; 


43 to dose stS3EELR> 62 yiseaesen esw  ti= 

‘+ lo doal .sutetaggs edd lo tnsbaegebal eboivs 

peel fat shies pnibesoxs 163 eefist Landtisaeqo as bad sco 
| er | pies 

mite® ¢lssioeaq o2 jud avisisqge sf 

fipeere 23 figaead saw oes \ tosmysteni doée 107? anoi2>dih 


~ a 


ie : 
aa 
a 


eyaqubens v1 erueesnd AG nol send (ee ty 


bi 
be 
ia jila> . Baw cesubengad--eweeerg fered / 
Lent Ssstic wisetdemscd . Noted) (sho Raa esnémurdekt fen xt 
tof Liiw One seemeizeni edt Asiw aidslisvs ei seu 3m 


i 
vivo | 6M idesdélao ens ewone f, + .e10g2% sted Osea oee 


ioad ont wootnem of Bean) +eoubenwad euteeerg.. sm 0.T 
tewllt $.39 seeerd siidw yn eae ods nidiiw guieaes 
gM Of ony esane . SeouBenerd Cl 0. S af} 203 ‘evi? a 
nese <stsmnoeiloznen ay oven ylioes iS JO1NsS ‘eoubshesd ia 
heidgere sd? 33 aes the Beitegs ar gino “aud oyfaes 1 


aa 5 
va kt ae 


7i sbealtapet ‘saw adijendhigs: “Agnotsibbe ne vsebedtis ai 


203 Ot stbeesng Tis isthaen ad gotaend2 Tag beviovat 
hire aoe) 

nogeiq seit no Se T19*9, bsof ay date vebnily> ais awstgets a 

=e! PAR 


BME. pNivong paabnede @ pike bemzotsed - ew ae 


| ; ; ee eel 
Mo (Pals ; bay hee 7 0 ; 
oa | ae eh gateat ai mwode “2 Ss 
oil 


Pg «lH +) Ni 


fog 


eG 


8e 


fie 


seonpsues, @iNSsaig edi OL JO} eAIND UOHeuqueg 1°D e4nbi4 


CAW) LAdLAG LINN NOTLISINO’ bLyd 


Oe 


9I ct 8 i 


LUaTIT4Ja09 UOTIeTaJJOI =J 
(2d) Ped’ =IdaosaluT A =q 
(AW/edW) 10ZZ"0 =auTT Jo adoTg =w 


(AW) indyno yTUN UOTITSTAbe e1eq =x 


(dW) aunssaud pattddy =) sayayy 
Q+Xl=\ 


(G666°0=4) SISA TUNY NOISSSY9SY YvSANT) 


(dd) SYNSSAYd Gal Idd 


ij 
\ 
i 
= \S) yy . " 
cs - 
=> ia : 
‘Sy es / 
; ¥ 
| 2 : , 
a 
} 
t Ls: 
he 
a 
= a 
- "= ~— 
a 
7 , - 
; ‘he 
\ 
! ‘ Sn 


ao | 


wr 


a 


GFP swan tallggh «¥ ers 


(Vel duarus jire naky a sta =x 
Veet? (65.0 -entl to sqale 
GRY MEOO~ wiqaowitnd Yd 


Pay 
, freind tee misvisren 1 7” 


deXan¥ 


ui 


el 


o"s 


o"e 


Onn 


190 


OF 


Jeonpsuesj SiNSSBid edi OT 40} SAIN uUOHeIquey z9 eunbi4 


CAW) LAdLNO LINN NOTLISINOY LUO 


g at 2- g- O1- 


JUSTIT 44909 UOTIeTaJJOD =J 

(Edy) 68ET'O =IdsovalUT A =q 

(AU/Ed) VZ6L°86-— =2UTT jo adoTs =0 

(AW) 1NdjNo 1TUN UOTITSTAbe e1eq =X 
(Fd) aunssaud pattddy =A sayayy 


q+XU=A 


(6666"0=4) SISATUNY NOISSSY9SY YWSNII 


ht- 


8t- 


Cc7 


o00¢c oost oooT 00S 


00Se 


(dd) SYNSSSYd Cal ddd 


ieee skaheearal : 


(iw sugrua simu raitlalupe ete x 
(Vee) ASSET BE~ “onli Wo eqole == j del 
: ‘DR BBEC0 ~tgaremt Y «d : 


< 


trsinatViens maltalerwo i : : 4 i 


ie 


O00T 


S28 


OS2 


JopuyjA any wesydeiq Jo} eAung uoHeiquey ED enbi4 
(be) YSONTIAD YIN OL SYNSSS3Yd GsaT Idd 


Se9 00S SCS 0S¢ Sel 


JUaTIT44Ja09 voTLeTaw0d =J 

(ND) €S00°0- = ydasuaiuT A =q 

(Edy/N>D 9SI0°O = BUTT yo adoTS =w 

(2d) JapuTTAS ste 0} aynssad pattddy =x 
(N™) JapuTTAD spe Aq paonposd peo =) taJayy 


8666'0=J = 4+X=A 


SISA TUNY NOISSSY9SY YYSNT) 


(ND) YSONITAD YI Ad da0ndOdd devo) 


‘ 


ee 


_ seniemtnine a — 


aa Sc sic te T 


: . é a eS ee Kein ani 


| | “B2O,0—7 dkny 
wa : 
GOO “abntlys 1 yd Geouberg bea) =¥ svar 
(G9) ~wonLlys als of srasserrg beligg? =X 
iDINIDD 8200 ~ gall No sqole 
' W0 SS8GO » seuretal Y J 
traialYtens neitaisym =1 . 
) | Ee | 7 
\ | ) wn 


38 


192 


Each of the above curves was plotted linearly and was 
subjected to a regression analysis to provide a good fit to 
the data points. The pertinent details of the analysis are 
contained within the figures. 

all linear regression analyses gave correlation 
coefficients of eed indicating excellent linear 
relationships. This implies data electronically collected by 
the data acquisition init can be converted with an extremely 


high level of confidence (negligible error). 


(ji) Calibration of Linear Varying Displacement Transducer 
mULVDT } 

The LVDT was calibrated with the use of a precision 
micrometer, the results of which are shown in Figure C.4. 
The regression analysis results are also contained within 
the figure. 

The linear regression analysis gave a correlation 
coefficient of 1.0 indicating an excellent linear 
relationship. This implies data electronically collected by 
the data acquisition unit can be converted with an extremely 


high level of confidence (negligible error). 


Calibration Apparatus 


The, apparatus. is «subject, to Significant Perror ~in 
measurement in four specific areas: vertical compression due 
to operation loading; vertical expansion due to operational 


heating; horizontal expansion due to operational heating; 
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and piston friction. To eliminate the influence of each of 
these factors from the experimental data a series of tests 


were performed to provide a set of correction curves. 


(i) Compression of Apparatus by Applied Loads for Different 
Temperatures 

A series of tests were performed to investigate the 
degree of compression the apparatus experiences during 
loading at different temperature levels. The specimen used 
in these experiments was composed of aluminum machined to 
the approximate dimensions of an actual oil sands sample. To 
eliminate any influence of induced pressure within the 
sample chamber, water was excluded from the test and the 
bottom drainage ports were kept open to atmospheric 
pressure. The temperature range for the testing was varied 
from room temperature to 200 °C and the consolidometer 
piston pressure from the 100 kPa seating pressure teaa 
maximum of 3500 kPa. 


The following test procedure was employed: 


e The 100 kPa seating pressure was applied at room 
temperature 

e The consolidometer was heated to the desired temperature 
level and allowed to stabilize 

e The LVDT reading was recorded 

® The consolidometer piston pressure was increased at 500 


kPa increments to 3500 kPa 


#e 


ce Tee 
os 


cy WwW 


* ~~ ° 
u , 
hid 2 
& aa & : 
7 
i 
> 
a ® 
J — Bi 
niet 
~ ote 
sale tor ie 
= = 
~~ © 


‘ 


iainemitegze eit mort atos26i sae 


Bi 
ss 
as 

W 


+tGe 


song Hesgeig uiestemobiloeno> efT © 


qoX si9w et36q ~speniath 


Ueeeig. ortigsse sti OOf > sft 


ot Beteset esw 1 


=e ves 


+arimifs. of eae 


e 


vo1g of bemrolisg ot 


rid ; mee 


yd euteraggh to nol eeenained Ww) 
aanutsneqng 


a eS 


i 


Oo petizs2 


cs ads nAoleestqmos ito 


s 
vA 
1utaecteo $nseisit:5 gs 
i o " F a aT" rie 26x85 
a - ~, * ¢. - 2s 
ad 2 enrenio.sJeamni xO190G0S 


on a ow 
: *Oo . so9naulinz Vie sisnimirt 7 
“vy a + « faToayw iedmado 


7 
s00S1 2iuIsisqmet ‘SAT 


1eqges MOOT Moy a5 
oof edz Ro. eiusesig » stig 
~N 
64% Q02E Bo uni xi 


giubscot¢ tas3 pniwolloi eiT i 


~ 


eivisisqmes 
etamobiloanes ait « 
ssilides2 o+ bswolls Bna fevel 
behtens2 2aw pnibses Tava ait # 


“Six 0 oae ost ainemsioni BTx 


195 


poe results derived! from! “this "seriesm@ o£. testing are 
presented in Figure C.5. 

To account for compression of the apparatus during 
actual testing the vertical displacement as measured by the 
tests at room temperature, 100 °C and 200 °C were averaged. 
The test performed at 50 °C was discarded due to lack of 
Similarity with respect to the other tests. From? Fagure: C,5 
zt is evident that influence of temperature on _ the 
compression of the apparatus by the applied loads was not a 


Significant factor. 


(ji) Vertical Expansion of Apparatus with Temperature 

The degree of vertical thermal expansion of the 
apparatus was determined for two seating pressures (100 kPa 
and 3500 kPa). As with the compression tests, the aluminum 
Specimen was used ana. the Sample chamber was held open to 
atmospheric pressure. The following test procedure was 


utilized: 


e The seating pressure was applied at room temperature and 
the LVDT reading recorded 

e Heat was applied to the consolidometer 

e After temperature stabilization occurred at the desired 
level, the LVDT reading was recorded 


e Steps 2 and 3 were mepeated yin 10) °C ancrements oy 200 2¢ 
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The pucees resulting from these tests are shown in 
Figure C.6 along with the theoretical vertical thermal 
expansion curve for the aluminum sample. The vertical 
thermal expansion of the apparatus is the difference between 
the expansion measured during the tests and the theoetical 
expansion of the aluminum sample. 

During experimental testing to account for the vertical 
expansion of apparatuS with temperature the results of the 
two tests were expressed by a linear relationship (Appendix 
D). Since the expansion of the aluminum sample is also 
linear, the vertical thermal expansion of the apparatus was 
taken as the difference between these two equations 


(Appendix D). 


(jii) Horizontal Expansion of Apparatus with Temperature 

The horizontal expansion of the sample chamber with 
temperature was taken into account by considering the sample 
chamber as an isotropic solid. When an isotropic solid is 
Subjected to a temperature change it expands with the 
distance between any two points increasing in the ratio a (a 
being the coefficient of linear expansion). This means the 
inside aianerer of the sample jacket enlarges in the same 
ratio as the external dimension. 

It can be proven theoretically that the coefficient of 
area expansion of an isotropic. solid) is "twree the 
coefficient of linear expansion (e.g., Halliday and Resnick 


.[1974]). Therefore, the change in area of the sample chamber 
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is 2 times the linear coefficient of linear expansion times 
the original area (at room temperature) times the change in 
temperature. 

The thermal volume change of the apparatus was taken 
into account in the measurements by multiplying the change 
in area of the sample chamber by the sample height. Details 


of the calculations involved are given in Appendix D. 


(iv) Piston Friction for Different Pressure Level and 
Different Temperatures 

The oedometer piston and sample jacket are composed of 
steel alloys differing slightly in material properties 
(Appendix B). As a ecul ee eee ee thermal expansion in 
the radial direction was of concern during operational use 
because of its direct effect on piston friction. Also of 
concern is the frictional resistance between the o ring and 
Sample chamber wall. It was therefore necessary to examine 
piston friction for both different pressure levels and 
different temperatures. 

The test procedure adopted involved holding the 
temperature constant and varying the pressure applied to the 
piston (Figure C.7). During the test the sample chamber was 
filled with water and the system totally deaired. The piston 
friction was taken as the difference between the pressure 
applied “to: thesepigpon (and the pressuremin, Cher watereas 


measured by the back pressure transducer. 
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The results (Figure C.7) indicate that the piston 
Eriction increases with both temperature and piston 
pressure. The maximum piston fraction pressure 
(approximately 230 kPa) was measured for the piston pressure 
of 3500 kPa at approximately 200 °C. 

Piston friction was accounted for in experimental 
testing by utilizing the piston friction versus temperature 


curve for the appropriate applied load. 
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APPENDIX D - DATA REDUCTION CALCULATIONS 


Volume of Extraneous Water 
Extraneous water refers to water not in sample chamber, 


le. water in porous stones, passages and ports. 


(i) Porous Stones 


TOD Top LOD. 
Stone (Inside) (Outside) (Outside) Bottom 
Weight wet St cw) 36.98 Silay ae) 95, 301 
(g) 
Weight dry 32.54 O2ece Sao 85.56 
(g) 
Weight of 4.66 4.69 Sli2o 9.75 
water (gq) 
Volume of 4.66 4.69 5125 9375 


wacer (cm) 


*Stone replaced after being damaged 


Total Volume of Extraneous Water in Porous Stones 


(a). Gnievally: 


4.66 + 4,69,49.75 = s'0 cm: 
(b) after top outside stone was replaced: 
4667+5, (25. 49.75. = 


19.53° cm? 
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(ji) Piston 


Thermocouple port: 


Thermocouple: length = 9.955 cm 


diameter = 0.318 cm 


Port: . length = 9.955 cm 

diameter = 0.397 cm 
Volume of extraneous water = Volume of port - Volume of 
thermocouple 


vome= 97 (9.955)'(0.25)1(0.397)* -— (0.378)? den0ed4s eme 


Port ccw of Thermocouple (viewing downward): 


length = 8.250. "em 
diameter ' = 0.318 cm 
length less plug length = 0.785 cm 
diameter of plug = 0.890 cm 


VOEk =era(0ii25 )d8e 250040508) trcigtGre2oimmess) (0. 890n2 =v e2 


eme 


Port cw of Thermocouple (viewing downward): 


length = 8.250 cm 
diameter = 0.318 cm 
length less plug length : = 0.860 cm 


diameter of plug = 0.890 cm 
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womet=) 7 (0.25) (82250) (05318)? = 2 (0.25) (8. 6a noTego) 2 


Total volume of extraneous water in piston 


= 0.443 + 1.142 + 1.188 = 


(iii) Pedestal 


Thermocouple ports: 


Thermocouple: 


Ports: 


Volume of extraneous water 


thermocouple 


Menem wm G7) 035) (Ones) (Oso 79% 


length 
diameter 
length 


diameter 


= Volume of port - 


(Oce48) se 90. 329 


204 


=1.188 


Ze Psocms 


S35 em 


316. cm 


LOS 57cm 


cme 


Back pressure drainage port including tubing to first 


Porto splug: 


WOR ean 0, 25) (3 2198) (O31 5) 


VOR = 700.25) (1.085)40 255) - 


length 
diameter 


Os Zoo wens 


length 
diameter 


Oey ecm 
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397 em 


Volume of 


valve: 
to5eacm 
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Plugs length = 1.435 cm 
diameter = 0.238 cm 


0.064 cm? 


VOE = 7(0.25) (1.435) (0.238)? 


Tubing: length = 8.905 cm 
diameter = 0.159 cm 


Weee=om(0.25)4(8.905) (0.759)? #90176 cm: 


Remaining drainage port: 
engin t=. 3.195)-om 
diameter = 0.318 cm 


G2 SeeMz 


VOE 


m(0.25) (3.195) (0.318) 


length less plug depth = 0.825 cm 
Giameter = 0.955 cm 


WOB- Sha (0225) 07825) 802.955) Reeder seapens 

Total volume of extraneous water in pedestal including 
tubing to first valve: 
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2.4630¢em" 


Total Volume of Extraneous Water; 
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Stone was Replaced 


Porous Stones (cm?) 19.10 19.535 
Piston (cm?) DIGI PAA 

Pedestal and Tubing 2.463 2.463 

(cm? ) 

Total (em? ) 24.336 24° 577.1 


Change in Sample due to Heating 


where: 

AA =change in sample area 

A = area of sample 

e® = thermal coefficient of linear expansion 
AT = change in temperature 


For the sample jacket (420 stainless steel) 


ae au10.8 x, 10° °/°C. (Appendix B) 


aere2 16 x 10° °7°C 


At room temperature (22 °C) the inside area of the sample 


chamber 18: 


Aveata( 764.2720" =" 4560.37 mm? 
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At any temperature (T,) the area (A,) is as follows: 
A; = A(1 + 2a(T, = 22]) 


Therefore, the volume (V,) at any temperature (T,) is as 


follows: 

WeeseA H- (1 + 2a[T, - 22]) 

Pre. evo ee4560naiahot iot epuératido ebT ,efu224 ml mmay 
Oe eV, =44560.37 Hit .0,0985 H(T, = 22) fmmed 


where: H = height of sample (corrected for vertical 


expansion with temperature) 


Therefore, the total volume of system at any temperature, 


aes 

Va. =< VORper4560037°H)+ 0.0985. (Ty -' 22) imma] 
where: 

VOE = volume of extraneous water 


T, = temperature 


H =height of sample 
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Vertical Expansion of Apparatus with Temperture 


From linear regression analysis performed on calibration 


data the following results were obtained: 


y =0.00448 x - 0.1324 


wheres 


vertical expansion of apparatus and aluminum sample (mm) 


<a 
i) 


temperature (°C 


“ 
u 


Considering expansion of aluminum sample: 


y' =0.00079 x’ - 0.0173 


where; 


y' = vertical expansion of aluminum sample (mm) 
x' = temperature (°C) 

therefore, 

een ay" 
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Yt = net expansion of apparatus due to temperature 


oeOn48)x--00.1328 390,000 79"x' 4 000173 


Ke 
u 


or; 
y. = OF0048ex6-n0.0007972k! 2=)0f4754 
but x = x', therefore, 


me = 0200201 xhe 0n1451 

Volume in Sample Chamber from LVDT Reading 

To calculate the volume in the sample chamber from the LVDT 
reading an aluminum sample was placed in the chamber at room 


temperature and at the desired ram pressure. 


From the LVDT calibration curve (Appendix C) 


where: 


y = data aquisition unit reading (volts) 
m = slope from linear regression analysis 


b = y intercept from linear regression analysis 


Since the height of the aluminum sample is known as well as 
the corresponding LVDT reading, the new y intercept can be 


calculated.- 


if eetonony ‘oa ¥ 2500040, - 3s 
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2.6 .@HD=1x26=P25525emm 

LVDT Reading for Aluminum Sample = y 

Therefore: 

b' = LVDT Reading - m (25.25) for aluminum sample 

As a result, the height of the sample during testing is: 


H = (LVDT Reading - b')/m 


Correcting for vertical expansion of apparatus with 


temperature: 


H =(LVDT Reading - b')/m - 0.0040 T + 0.1151 


and as before: 


Maw= VOE#+4650537 H+ 0.0965 H (TT 9 2220) 


Calculation of Applied Ram Pressure 


From calibration curve for 2.0 MPa transducer: 


y =mx +b 


where: 


we 
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y = applied pressure (kPa) 

x = data aquisition unit output (mv) 
m = slope of line (kPa/mv) 

b = y intercept (kPa) 


Diaphram air cylinder: 


wheres: 


ve load produced by air cylinder (kN) 


x' = applied pressure to air cylinder (kPa) 
m' = slope of line (kN/kPa) 
b' = y intercept (kN) 


Ve 


Since y = x 


vane ane MCX ct Dae, 


Therefore: 


The applied room pressure = y'H(1000)7/V, + Weight of piston 


and concrete cylinder (kPa) 


where: 
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22 
V; = volume of sample chamber corrected for extraneous water 


and radial expansion of sample jacket 


H = height of sample corrected for vertical expansion of 


apparatus with temperature 
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APPENDIX E - COMPUTER PROGRAMS 


2 KK ok oR ok ro Kok oF ok kK ok ok ok ok ok ok ok oko KK ok ok ok ok ok kok 


PEMD2 eT 
VERSION 1.1, MAY 1972 
TWO-DIMENSIONAL HEAT CONDUCTION 
USING 


FINITE ELEMENT ANALYSIS 


* * 
* * 
* * 
* * 
* * 
* * 
* * 
* * 
* * 
* * 
* * 
* * 


2K 2K KK OK OK KK oO KK oko ok OK ok 2K ok ok ok ok ok ok ok ok ok oF 


eK OK KK RK KK oo ok ok ok ok ok ook ok ok ok Oo Ko ok ok ok OK ok ook eK oR ok ok ok ok ok ok ok ok Ok ROR EK OK Kh KKK KEKE 


FEHT2 (FINITE ELEMENT HEAT TRANSFER-2 DIMENSIONS) 1S A PROGRAM 
FOR THE ANALYSIS OF TWO-DIMENSIONAL HEAT CONDUCTION PROBLEMS BY A 
FINITE ELEMENT TECHNIQUE USING LINEAR TRIANGULAR ELEMENTS. 


THE PROGRAM WILL HANDLE THE FOLLOWING PROBLEMS. 
1) STEADY STATE PROBLEMS. 
2) TRANSIENT PROBLEMS BY THE FOLLOWING METHODS. 
A) EULERS METHOD. 
B) THE CRANK-NICHOLSON METHOD. 
C) A PURE-IMPLICIT METHOD. 


THE FOLLOWING BOUNDARY CONDITIONS CAN BE CONSIDERED. 
1) TIME INDEPENDENT CONVECTION. 
2) TIME INDEPENDENT HEAT FLUX. 
3) TIME INDEPENDENT TEMPERATURE. 


THE FOLLOWING PROPERTIES MAY BE GEOMETRICALLY VARIED. 
1) THERMAL CONDUCTIVITY. 
2) DENSITY. 
3) SPECIFIC HEAT. 
4) RATE OF HEAT GENERATION PER UNIT VOLUME. 


IT IS ASSUMED THAT THE MATERIAL IS ISOTROPIC. ANISOTROPY CAN BY 
CONSIDERED BY MAKING AN APPROPRIATE GEOMETRIC CHANGE OF SCALE. 


THE PROGRAM IS DIMENSIONED TO ACCOMMODATE SOO NODAL POINTS 
PROVIDED THAT THEY ARE NUMBERED TO RESULT IN A BAND WIDTH OF 30 
OR LESS. THE BAND WIDTH IS EQUAL TO 1 PLUS TWICE THE MAXIMUM 
DIFFERENCE BETWEEN NODAL POINT NUMBERS OF AN ELEMENT. 


ANY CONSISTENT SET OF UNITS MAY BE USED. 


REFERENCE- 
MYERS, G.E., (1971), ANALYTICAL METHODS IN CONDUCTION HEAT 
TRANSFER, MCGRAW-HILL BOOK COMPANY, NEW YORK, N.Y. 
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Fe OR OK KK KOK OK KOK KOK RRR OK ORK OR KOR ok kok ok oF ok ok ok ok ok ok ok ok 


DIMENSIONS 

ROR ORR ROR IRI GI aR kk ak dca ak ae ak ok ae ak 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON BN(500, 30) ,CN(500, 30) ,SN(500,30),SP(500),Y(500),X(500), 
*T (500) ,GN(500) ,RN(500), GEN(50) , CEE(50).RHO(50) .CAY(50),TIMAX, 
*DTIME,TIME,NTOUT (500) ,MN(500) ,NRL(500),KN(500),JUN(500),IN(500),uJC, 
* ICOUNT ,NTIME,NT,NQ,NH,NC,NG,NE,NN,IBW,IOUT,IIN 


Dk KK kK oo ok OK ok oo Ko OK ok 2K ok ok ok ok oR ok Kok oko oR KK KK a KK ok ok ok ok Kok ok ok ok ok OK kok oR ob ok ok ok ok ok ok 


SET INPUT/OUTPUT FILE NUMBERS AND INITIALIZE 

EEK ROR OR OR ROR IO RR OOK ER OR IG RoR aR a ok Ga ak ak a a ak ak ac af ak 
LUNG ==5 

IOUT = 6 

LABEL=O 

IFLAG = O 


DK OK EK KK KOK KK KK KK OK OK KK ok kK ok ok KOK OK ok ok OK OK ok ok Oo OK OK OK Ok oF ok ok ok KK ok OK OR Oo KOK OK KOK OKO OK Ok Ok 


MAIN PROGRAM 
SKK KK KK KK ok ok KK KK Kk ok oR ko OK ok Ko OK OK ok OK 2k Oo ok ok ok ok oR ok oR OR EK oR ok Oo OK oo OK ok ok kk ok KOE OK OK Ok OF OK 
CALL DATA(LABEL) 
IF (LABEL.EQ.1) GO TO. 15 
CALL FORM 
PEMENE) 02, 273 

IF (NQ) 4,4,3 

CALL BINT 

IF (NC) 6,6,5 

CALL INIT 

CALL TRANS 
TESONT) $8,807 

CALL MODSN 
CALL DBAND(NN,IBW,SN,LABEL) 
IF (LABEL.EQ.1) GO TO 14 
CALL RHS 
PEDENT) 140474210 

CALL MODRN 
CALL SBAND(NN,IBW,SN,RN,T) 
CALL RESUL(IFLAG) 
KE LONG) tae t ti2 

IF (TIME - TIMAX) 13,1,1 
IF (IFLAG) 16,9,16 

WRITE (IOUT,601) 
STOP 
WRITE (IOUT,602) 
GO TO 15 


Se KOK ROK OK OR ROR RO RO OK oR EK ok ok oko oR oR ER RR Bo oR Rk oa oR oR Ro oR oR oR ook ok Bok oe oe ok 


FORMAT STATEMENTS 


eK OK OK Ko ok ok ok ok ook ok ok kok ok ok ok ook ko kK oo ok Ko oo ok ok ok ko ok kK ok OR OK OK OK OK OF 


FORMAT (//’0’,15X, ’DECOMPOSITION FAILS’//) 
FORMAT (’0’,15X, ‘OUTPUT TIME SPECIFICATION ERROR’ ) 


END 


214 


ss 


oe bios elses es eegede * 


ie 


enatnae Poe eneend es Peet ee n° 


tobe tqoees, (oa : Bi tate liter ik ky ASE 70 + 1 —s =™* 
SAME? 408) baa, £08 MINORS 139. (02 MBO (008 J et 
2, (RAE NL. 1008 IML (062 SL (G08 )JaM. (OGr Wel, (0o8} TG 
Wit, TUS War 


“TID tN e he eee eb Oe Tee ee he eee Oe 


a 


CED OFOl OOOO CAO Or Oi@i@: OI@1® 


Or@ 


N9NNdNaN0= Aoa ao CY AKO OF Qi OE AO 


Qn 


SUBROUTINE DATA(LABEL) 
IMPLICIT REAL*8(A-H,0-Z) 


eK OK KK OK OK KK ok ok ok ok 


* READ INPUT DATA * 


We KF KKK KK OK KK OK OK OK OK OK 


FEOF RE KK KK KK KK OK RK KK oR KOK OK KK OK OK RK OR OR KOK OK OK oo ok Kk ok ok ok ok ook Kok ok ook ok ok ok ok 


ARGUMENT - 


EABEE SEND OF VEIL EES INDICATOR 
="© IF BXECUTION IS TO CONTINUE 
= 1 IF EXECUTION IS TO TERMINATE 


2K oe ok KKK OK ok ok OK OK oO KK KR oo Ko ok ak ok ok ok ok kok ok KK oR KO oR oR ke ok ok ok Ok 


De KK Ro ok OK Ko ok KK oo OK KK oo OK ook ok ok ok ok ok ok ok ok ok dk a ok oo oo ook ok ck of ok ok ok ok ok 


DIMENSIONS 

KK KK OK KK ok ok OK oF KK ok ok ok OK Kk kK OK Kk OK OK OK OK KK ok KO oF ok OK OF OK OO KOK OK OR OK OK OK KOK KKK 
COMMON BN(500,30),CN(500,30),SN(500,30),.SP(500),Y(500),X(500), 
*T(500),GN(500),RN(500), GEN(50),CEE(50),RHO(50),.CAY(50),TIMAX, 
*DTIME, TIME,NTOUT(500O) ,MN(500) ,NRL(500),KN(500),JUN(500),IN(500),uUC. 
* ICOUNT ,NTIME,NT,NQ,NH,NC,NG,NE.NN,IBW,IOUT,IIN 


DIMENSION TITLE(20) 


eee ee ee ee ee ee ee a ee ee ae a ee ae ae ae ek oe ee ae Ce aed 


PROBLEM IDENTIFICATION 

KOK KK OK OK OK OO oo RO OR Fo ROKK OO oF OO OK KOR ROR Kk KOR OK KO 
READ (IIN,501,END=99) (TITLE(I),1=1,20) 

WRETE (LOUTS COM) CLUMEE GI). 1=45,20) 


we KR OK KR kK ok ok ok ok OK KOK ok OK KK KF OK Ok ok ok ok kK eK kK ok ok OK ok KOK Ok ok KK OK ok ok ok ok ok OO 


PROBLEM INDICATORS 

sk ok ok ok Kk ok oe ok OK OK KOK KK OK KK ok EK OK OK OK OO OO KO Ok OO kk oe 
WRITE (IOUT,602) 

READ (IIN,502) NG.NC,NH,NQ,NT.NN,NE,NM 

WRITE (IOUT,603) NG,NC,NH,NQ,NT,NN,NE,NM 


ORK OK OK OK oo oo Ro oR Ko oO oO RK OK OK oR OR oR oR bo OK OR Fok oh ok ok oO 


MATERIAL PROPERTIES 
OO KK KK KF KK OK ok ok ok oR ok OR OO RO OY OK Fo OK Ok oh OF 
WRITE (IOUT,604) 
DO 1 K=1,NM 

READ (1IN,503) M,CAY(M),RHO(M),CEE(M),GEN(M) 

WRITE (IOUT,605) M,CAY(M),RHO(M),CEE(M),GEN(M) 


Seo OK KO OR Ko oo ko ok oo oe eo OK oo oo OK oR OK ok oo OF oR oh RR oP ok oO 


NODAL COORDINATES 
Se EK OR OK OK RK RK ook oo do oF oR OR KOK OR KOR OO ROR KK EE OO 
WRITE (IOUT,606) 
DO 4 K=1,NN 
READ (IIN,504) J,X(J),Y (uv) 
WRITE (I10UT,607) J,X(JU),Y(U) 
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HEE KK OK ko OK KK KK Ko OK ok KOK OK KK ok KK Rk kk ok ok KR KOK OK KO OK KKK Kk kk ok ok ok ok ok ok dk ok 


ELEMENT INFORMATION 
2K KK KK KK KK ok OK ok ok oR KR ok oR KK KK KKK ook ok ok oko ok kook RK ok ok ok ok ok kok kok kok 
WRITE (IO0UT,608) 
IBW = O 
DO 11 K=1,NE 
READ (IIN,505) I,IN(I),UN(I),KN(I).MN(I) 
IF (IN(I)-UN(I)) 6,6,5 
ITEM = IN(T) 
IN(T) JN(1) 
JN( I) ITEM 
IF (JN(I)-KN(I)) 9,9,7 
ITEM = JN(I) 
JN( 1) = KN(I) 
KN(I) = ITEM 
IF (IN(I)-UN(I)) 9,9,8 
ITEM = IN(I) 
IN(I) = JN(I) 
UN(I) = ITEM 
ICMAX = KN(I) - IN(I) + 1 
IF (ICMAX-IBW) 11,11,10 
IBW = ICMAX 
WRITE (IOUT,608) I,IN(I).UN(I),KN(I),MN(I) 


PABEE = 70 
RETURN 
BAB bie sa 
RETURN 


OK KK OK OK OK OK OK KK KO KK KK Oo OK OF RK ok Ok KK OK OOK OOK OK OK OK OR KOK KO KO Ok OK 


FORMAT STATEMENTS 

OK oO OK OK ok ok Kk KOK OK ok OK OF OK OK KO OK OOF OF OO OR OK KOK OK OOK OK OF OK KOR OK KO OK OK OE OF OK KO KORY OO Ok OF OF OK KOK 
FORMAT (20A4) 

FORMAT (815) 

FORMAT (15,4F10.0) 

FORMAT (15,2F10.0) 

FORMAT (515) 

FORMAT (/1’,15X,20A4) 

FORMAT (’0’,15X, ‘NG NC NH NQ NT NN NE NM’ ) 


FORMAT (”% ”, 15X,12,915) 

FORMAT (’0’,15X, ‘MATERIAL PROPERTIES’/ 

1 TEAS Oma ATi le: CAY RHO CEE GEN’ ) 
FORMATIN( lee fox. 4h a2) 

FORMAT (’0’,15X, ‘COORDINATES OF NODES’/ 

{ OU TAGE ONODE 0 4K. 7X? Bx Vs ) 

FORMAT (° °’,415X,13,2F9.3) 

FORMAT (’0’,15X, ELEMENT INFORMATION’ / 


1 Re NO Xe ae At. I J KO MATIC 4) 


FORMAT 1( 4 © 715% F139 16 254, FS) 
END 
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SUBROUTINE FORM 

IMPLICIT REAL*8(A-H,0-Z) 

OF eK oo oo oR KK ok KK ok ok Kok ok ok ok ok kok ok ok ok ok ok ok kok ok 
* FORM ELEMENT CONDUCTION, GENERATION, * 
* AND CAPACITANCE MATRICES * 


Ae EK OK KK KK KKK KK ok KK OK OK ok ok Kk ok ok ok KK ok ok ok OK 


2K KKK Ko KK KR KK KOK OK KK KK oo oe ok ok KR KR Rok ok ok oF ook ok ok ok ok ok OK Ko a EE kK ok ok ok 


DIMENSIONS 

3K KOK Ke ok ok ok kk KK ok ok OK ok ok Kar OK ok oo ok OK ok ok ok ok ok ok ok OK ek ok ok Ko ok kok kk kok ok ok kok ck ok ok ok kk ok ok ok 
COMMON BN(500,30),CN(500,30),SN(500,30),SP(500),Y(500),x(500), 
*T(500) ,GN(500) ,RN(500) ,GEN(50) ,CEE(50),RHO(50).CAY(5O),TIMAX, 
*DTIME , TIME ,NTOUT(500) ,MN(500) ,NRL(500) ,KN(500),UN(500), IN(500),uC, 
*  ICOUNT,NTIME,NT,NQ,NH,NC,NG,NE,NN,IBW, JOUT,IIN 


DIMENSION ID(3), SE(3,3) 


eK OR OK KK KK OK KOK ok oF ok KK OK OK ok Ok KK OK kK Ro ok KK Oo OF OK OK KE OK OK OK OK OF KO OK 


INITIALIZE 


OK OK OK OK OK OK oF OK KKK ok Ko KOK KK OOK KO RK OK KK ok OK oo ok oF ok OK KK OK OK Ok ok ok Kok ok KF OOK OK ok oF ok ok 


DO 1 J=1,NN 


eK OF KK KK OK OK ok ok oR KK OK OK OK oe oR Ko EK of KK oo Kk Ok ok OK KOK OK OK OK OF OK Ok Ok OK Oo Ok 


ELEMENT CONDUCTION MATRICES 
OK ok Oo Kok OK ok OK ok ok ok on oR OK OK OF OK KO OK ook OK OK OK OK OK ok OK OOK Ook OF OK OK OK OOK KO OK OY OOF OY KOK KOO OK OK KO KEE 
I fe) 
I Loctet 
IF (I-NE) 5,5,99 
INN=IN(I) 
JNN=JUN(T) 
KNN=KN(I1) 
MNN=MN(1) 
XI JU=X(UNN)-X (INN) 
XIK=X(KNN)-X (INN) 
X JK =X (KNN)-X (JUNN) 
YIJU=Y(UNN)-Y (INN) 
YIK=Y(KNN)-Y (INN) 
YUK=Y (KNN)-Y(JNN) 
AIJK = DABS(XIU*YIK - XIK*YIUJU)/2.0 
E2=(CAY(MNN) )/(4.0*AIJK) 


SE) See GXUK 227 VK 2p) 
SE(1,.2) = sE2*(X1IK*XdK + VIK*VYUK) 
SE( 1,39) = E2*(XIUtXKUK + VIU*VUK) 
SEC(202)) = EltCKIK ts? ee vIKer2) 
SE( 23) = -E2@( Xi UEKI K FOV TdF VIK) 


SEG 53). = EQFCKItt2 + ViIdt*2) 
ID( 7) = INCI) 


TH(2) = INCE) 
ID(3) = KN(I) 
DO 6) J=1,9 

TR: = ID(J) 
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DO 6 K=J,3 
TCe="1DiCK)2 — REA 
SN(IR,IC) = SN(IR,IC) + SE(JU,K) 


ROK OK KK Ko oR OK OK ok OK OK oo KK KK oo ok ok oF ok ok KOK oR ok OF KK ok oo ok oo Ok ok OK OK OK KOK OK OK OF OOK OK OK OK OK KO 


ELEMENT GENERATION MATRICES 
2K > KOK kK OK OK OK 2K OK ok ok KO OK OK OK Ko KOK OK OK OK OK Kk OOK OK OK OK OK KOK OK OOK OK OK Ok OK OK OK OK OK OK OO OK OK OK OK OK OK OK Ok OK KOK OR OK OK OK Ok 
IF (NG) 9,9,7 
E2=(GEN(MNN)*AIUK)/3.0 
DO 8 J=1,3 
IR’ = ID( Uv) 
GN(IR) = GN(IR) + E2 


DKK KK OK ok oR ok ok ok ok ok ok ok ok KK kK ok ok ok ok Ko Kk ok ok ok ok ok ok ok oh ok ok ok ok ok OK KOK ok Kok ok ok ok ok ok OK oR OK OK ok ok ok ok 


ELEMENT CAPACITANCE MATRICES 
SK OK OK KK KK ok Ko Ko oF OK OK OK KK OK OK oh OO OK OK OK OK OK KOK OK OK OK OO OK OK OO KK OK RK OK OK OK kK OK Ok Ok 
1S. (ONG). W265 122, VO 
E2=((RHO(MNN) )*(CEE(MNN) )*AIUK)/12.0 
SECM, 1/)' @ 12 sOFE2 
SE (A 320i =4E2 


SE(423) =) £2 
SE(Q, 2 2 .OFE2 
SE (2, 3)' =" £2 
SECS; S08 =h 200*E2 
DO 11 d=, 3 

2 ss. say(ul) 


DO AdeK=U.3 
Ic = 0b (K) ei wR) + 1 
CN GLRATC )- =F ICNCER: 1G) cee SE CuK) 
GO TO 4 
RETURN 
END 
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SUBROUTINE BINT 
IMPLICIT REAL*8(A-H,0-Z) 

DEK ER KK KKK Ko KOK OK ok ok OK OK ok ok ok ok ok ok ok KK ok ok ok ok ok ok KF ok ok ok ok Ok xk 
* SET UP MATRICES FOR CONVECTION AND NON-ZERO * 
* HEAT FLUX BOUNDARY CONDITIONS x 


2K OK OK KK ko KK KKK OK KK oo KOK ok ok KOK KK KOK OK ok ok ok OK OK ok ok ok OK ok ok Ok 


2K KK KE KK Oo KK KKK OK oF KK KK ok ok Ko ok ok ook ok oo KK ok ok ok ok ok ok ok ok ok ok ok ok ok ok 


DIMENSIONS 

ROK OK OK KK KK KOK KOK OK OF OOK OK OK KK OK ok OK OK KK KOK oo ok ok ok ok ok ob ok ok ok oF ok ok ok KOR KO KOR OK KOKO EK KOK KARE 
COMMON BN(500,30),CN(500,30),SN(500, 30) ,SP(500) ,.Y(500),X(500), 
*T(500) ,GN(500),RN(500),GEN(50) ,CEE(50) ,.RHO(50),CAY(50),TIMAX, 
*DTIME,TIME,NTOUT(500) ,MN(500) ,NRL(500) ,KN(500),JUN(500),IN(500),uUC, 
* ICOUNT ,NTIME,NT,.NQ,NH,NC,NG,NE,NN,IBW,IOUT,IIN 


We KK KK OK KK KK KK KKK KK kK ok Ko ok ok kK ok oF Xe ok oo ook ok KK ok ok ok ok ok ok ok ok ok Kok ok OK OK OF OOK Kk ok 


CONVECTION BOUNDARY LINES 
SEO KOK OR ROKR Eo EO oR Rok ok KR ok OK RK Kk R OR Ok orb kk 
ile (NED) S54 S)5 4 
WRITE (IOUT,603) 
DO 4 J=1,NH 
READ (IIN,502) IS,US,H, TAMB 
iN? (NS = US) S82 


TTLEMe= LS 
IS = JS 
US ee 


WRITE (I0UT,604) J,1S,JUS,H,TAMB 
NG US GES 


XIU = COS) =) KCTS) 
Wiig) Ss Ws), = AEs) 
SIU = DSQRT(XIU**2 + YIU**2) 


E2)=-(H*SIY)/6.0 

SNGES th = SN(TS, 1)" 2.0782 
SN(IS IC). = SNOIS.ICY) + (62 
SN GUS Hal): = sSNGUS al st 2 ORE 2 
E2 = (H*SIJU*TAMB)/2 


GN(IS) = GN(IS) + E2 
GN(US) = GN(UJUS) + E2 
CONT INUE 


OK OR oR OR Kk ok ok ok Kk ok Ro KK OK OO KK KK OK ok OR KOR OK KOK RE oo KK Ko OK OK oR OR OF OK OR KO OR KOK KO OK OK 


NON-ZERO HEAT FLUX BOUNDARY LINES 
eK OK OK KK KK KK OK oR Ko ok ok oo ok ok ok KK oR ok OK oo ok oF RK Oo OK ok ok Ek Kk Kk ok ok ok 
IF (NQ) 10,10.6 
WRITE (IOUT,605) 
DO 9 J=1,NQ 
READ (IIN,503) IS,uJS,Q 
EF CUS’ = JS)" 8 ee 


ITEM, = IS 
kSe=euUs 
JSS SS WSh 


WRITE (IOUT,606) J,IS.US,Q 
Sai) a SCS). = 36S), 

Witul = Wiis) = YCus) 

SIU = DSORT(XIJU**2 + YIuU**2) 
E2 = (Q*SIuJ)/2.0 

GN(IS)> = GN@iS) + 7&2 
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GN(JUS) = GN(JS) + E2 
CONTINUE 
RETURN 


OK OK eK KK OK ok KKK ok OK OK OK OK OK KOK OK KKK KK KK KK OK KK OK KK KK KK OK ok Oo OK Kk oF Ok ok Ok OK Ok 


FORMAT STATEMENTS 

KK KK KK KO OK KK KK oR OK ok KOK ok ok ok ok ok ok KOK ok OK OK ok ok ok OK ok OK OK OK OK OOK OK OK OK OK OK OK OK KOK OK OK KOK KK OK OF OK OK OK OK OK OK 
FORMAT (215,2F10.0) 

FORMAT (215,F10.0) 

FORMAT (’0’,15X, “CONVECTION BOUNDARY LINES’/ 

4 ro 45X, "CINE I J HCON TAMB’ ) 
FORMAT (’ ’,15X,13.16,14,F8.2,F10.3) 

FORMAT (’0’,15X, ‘NON-ZERO HEAT FLUX BOUNDARY LINES’/ 
{ 1 EAS AC IINE I rv) FLUX) 

FORMAT (’ ’,15X,13,16,14,F8.2) 

END 
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SUBROUTINE INIT 
IMPLICIT REAL*8(A-H,0-Z) 


eK RK RK KK KK ok ok ok ok OK Ok ok ok ok Ok 


* READ INITIAL CONDITIONS * 


EK RE KK OK OK KO KOK KOK KOK OK OK OK Ok OK KOK OK 


DEO KK KK OK OK OK OK KO KK KK OK Ko KK OK KK OK FO Ok ok kK KK OK KO oF OK OK OK ok OK OK KOK OK OF KOK 


DIMENSIONS 

eK KE KK KE KEK OK EK KK KK EK ok KK KK ok ok ok ok oe oR ok ok ok ok ok OK OK OK ok OK KK OK OK ook 
COMMON BN(500,30),CN(500,30),SN(500,30),SP(500).Y(500),X(500), 
*T(500),GN(500),RN( 500) ,GEN(50) ,CEE(50),RHO(50),CAY(50),TIMAX, 
*DTIME, TIME ,NTOUT(500) ,MN(500) ,NRL( 500) ,KN(500),UN(500),IN(500),uUC, 
ICOUNT ,NTIME,NT,NQ,NH,NC,NG,NE,NN, IBW, IOUT, IIN 


WOK KK Ko ok KK KOK KO KK oR OR KO KOK OK oo KK ROE OK OK ok ok ok ok ok ok ok ok ok ok ok ok ok ok OK KK OF Ok OK 


INITIAL AND MAXIMUM TIME 

SE OK KR RRR OR EE OK OK OE ORE OR RO dR OR OE kk oo oR A Rak oR deat ok ok 
READ (IIN,501) TIME,TIMAX 

WRITE (IOUT,601) TIME, TIMAX 


DK KK KK Ro KE KK oo KR OK OO oR OR OK OK OK OK OK OK OK OO KK KOR YOK OK OK OK OK OK OK ok OO OK OK ok Ok OF 


INITIAL CONDITIONS 
OK ok ok ok ok ok kK KO OK OR OK OK OKO OK KOK oR KK OO OK OK KK KK KO OO KOK OK OK KO OO OK OO OK OE Ok ok 
DO 1 J=1,NN 
READ (IIN,502) K,T(K) 
WRITE (IOUT,S02) K,T(K) 
RETURN 


OF OO OF OK ok ok ok ok ok ok ok ook KK OK Ko OK OK Oo KK OK KK OK KO KO ok ok ok tO OK ok ok OK OF ok Ok Ok Ok ok 


FORMAT STATEMENTS 
sk ok ok ok oR OK KR KK KK KK OK OO OK OOK OK KO OK OK Kk OK OK OK Ok KOKO OK KK OK OO OF OK OK OK 
FORMAT (2F10.0) 
FORMAT (15,F10.0) 
FORMAT (’0’,15X, INITIAL TIME=’,E£13.6,4X, “MAXIMUM TIME=',E13.6/ 
1 'O’,15X,’/ INITIAL TEMPERATURES’ / 
* +, 45X, ‘NODE TEMP.’ ) 


FORMAT ( 9545xX, fa 3FH0. 3) 
END 


22H 


NMMNNMNRHTTINE Re (2 cr rsaSe 


Ses ot ween nn eles eee 
da 


108%, (0027 Y . (ORR Tre 
SREY 1687 TAS (OE) 


4 a 


Oe hes 


Sia J 
ou, (S08? . tone ae pt 


“ii 


il “a 
1 


ois 4 0 04 Ce meme si 5 


€*+ e@n Ce 


eee =< @ o's 


Sei ene oe 


vas oJ 


\ sts Rweees oped - tetowewt orton aes; 


9 erg, 


er 


A I — fr al 
me o4e4 eset odes & Feber tr ee ASPEN Es eae a ewe 


; “ary 
= 26: ¢543 Op bee & @ ni aerate seers t= 


Naat aria ,. 


Pa ¢ al 


Hivensrinewartevinheta sda at o- eee 
nel 
lth na 


, T 
yt 
Aa tH tees ar Ho eeean tame unle enensinad pave 


“RIT winks ea 9S," 


iG Ore) rar Gra 


HOLOare 


A QGiOral@ 


OLOr@ Ola 


De) =AAQQ0N90 


WAIANQAAA 


SUBROUTINE TRANS 
IMPLICIT REAL*8(A-H,0-Z) 


EK KKK KKK KKK ok oR OK ok oF ok ok kK ok Ek ok ok ok kk 


* SET UP TRANSIENT ALGORITHMS * 


DK OE KK KK KK KOK OK OK OK OK OK KK oF OK KOK KOK OK 


FEO KK OK ORK KK KR ER KR ROKK KO GK KOK KR aK Kok ok ook kok 


DIMENSIONS 

eK KOK OK ok ok Kok KK OR Kok KK ok a KK ok KOK a oR ok ok ok ok Kk oo ok ak ok ok oR OK oR oR ok ok ok ok ok ok ok ok ok ok ok ok ook 
COMMON BN(500, 30) ,CN(500, 30) ,SN(500,30),SP(500),Y(500),x(500), 
*T(500),GN(500) ,RN(500), GEN(50) ,CEE(50),RHO(50),CAY(50),TIMAX, 
*DTIME, TIME ,NTOUT (500) ,MN(500),NRL(500) ,KN(500),UN(500).IN(500),uUC, 
* ICOUNT ,NTIME,NT,NQ,NH.NC,NG,NE,NN,IBW,IOUT,IIN 


OK KKK OK OK KOK KK KK KKK KK KK oO ok ok ok ok ok ok ok ok ok Ko ok ok OR OK ok Ok oR ok OK ok Ko ok ok ok oF ok ok 


TIME INCREMENT 

KK KK OK KK KK KK ok eK oo RK oo KOK Ko ok oe KK OK KK OK OK OK Ok KK KO Ok oo kk ok KO ROK OK KR OK OR OR OK OK KOK 
READ (IIN,501) DTIME 

WRITE (IOUT,601) DTIME 


WK KK OF Ko KK KK OK OK KOK OK OK KK OK kK ok ok kK ok ok ok ok ok Ok kK OK OK OK OK ok ok OK OK Ook OK OR ok KK OK KO OOK OF KO OK OOK 


PRINTOUT TIMES 

OK OK KK OK KO KOK KOK KK KK KF KK OK EK OR OR OO KK KO OK KK KK OK OK OK OK OK 
READ (IIN,502) NTIME 

WRITE (I0UT,602) NTIME 

READ (IIN,503) (NTOUT(1I),1I=1,NTIME) 

WRITE (IOUT,603) 

WRITE (IOUT,604) (NTOUT(I).I=1,NTIME) 


ak eK KO Kk KK ok ok ok ok OK ok ok ok KK ok ok kK ok OR Kk ok KOK oo OK oo ok ok Kk ok ok ok ok ok ok ok ok OK KO OK OK KOK KOK Ok Ok 


TRANSIENT COMPUTATIONAL METHOD 

Se OE ER OR OR Ro BR ooo OK OR OO OR RR Rk ok ok kor 
ICOUNT = O 

UG =o4 

GOMNOm GiGi on eae 


OK OK KO KO OR OK ORO OK OK OK OO KO EK KR ok OK KK OK oR KR OK OR OR OR OO KO RO OK KO OK 


EULERS METHOD - NC=1 
OK oR OK KK kK oo Ro KK oR KOK EK KK OE OO KK KE KO OF OO OK OK KOK FO 
DO 2 J=1,NN 
DO 2 K=1,IBW 
BN(J,K) = CN(JU.K) - DTIME*SN(U,K) 
SN(U,K) = CN(JU,K) 
WRITE (IOUT,605) 
GO TO 7 


Se EO KOR ER OR I RRR RKO a KR KK RK ro ok 


CRANK-NICHOLSON METHOD - NC=2 
eR OK EE ROR OO eo OR I RIOR OE ROE oR oR oR oo EE EO OF ok oe kok ok ok ok ok ok ok 
Diti2Q= Dil ME/ 2-0 
DO 4 J=1,NN 
DO 4 K=1,IBW 
BN(JU.K) = CN(U,K) - DT2*SN(U.K) 
SN(JU,K) = CN(JU,K) + DT2*SN(J,K) 


Ze 


oi y ; 
a a) : 
} a i" i) 
+ ree ’ . 
bs te f 2 ce 
: 7 <i onal eae 
it re en ae | 
 APTRET,, £Ge ro) OLE Monae (ahaa . Wagers: 
G1 io ot (od. Ka, ese 
“tl. Tint we sia eee setiall oh, 
3 vr oa 
; |_| aa eiersdnaaen ihe hmgial + knt hb GS serevcsene , 


‘envwete sco seewhene oT : 

“gaia Poe wit) daag 

i) Ms 6 ett nen FOprT - mM 
: ' 4 


rae seek =. MaMasontr- psi 


Pe or “ae wreaks Lenehan signe au othe 

Abe JIT TUOTHI9 

ack Ys medgers yedyae a hela :< carayateegh ee aousesesenes, = 

. : ef aMaTy (Boa Wit) Gag R= ar 
ree re 

{wT ott. “(yeu 6a ms Po fd 


ese ae a 


aie: Vets lipesset asia 


* 
” ‘ ‘ a 1 A 
bawe &4 90h 6 oS Vite ea, CA eae or aaa 
8 inate pac a 


7 A= acer ee ee es a nee 44+ 98 69 F eee Ta, Sea 
\ . s 7 - 2 ( : 5 - A his oe 
a: a x) ie: a thie 3 m3 te 

; a : a) ve vf es ; 
or 2 : > a “fd fe 


; 4 “¥ Z 
oe bo we tld Havisceintiity he 
7 . | as 
T 


? 


A 
fn eeeren ees ‘h4 7 : * : i Hoes ss = 
a ices ae | OBHTS | MOORS Taha 
: 7 aweee se Rea > re : “ty ia oir nk ae ; 
7 ; : : : . a Fs m ors : ea 
i 7 , / - * 


aandgaand 


o 


WRITE (IOUT,606) 
GO iG) 87 


OK OK KK KK KK Ko oo OK KK ok ok ok ok ok ok ok ok ook OK ook ok ek ok ok ok ok ok ok ok ok ook ok ok ok ok ok ok ok ok ok oF kok kok ok ok kok 


PURE-IMPLICIT METHOD - NC=3 
EK OK RK oo Eo OK ESO ORE Rk ok ko Sak ao of ak kak ok dee tek 
DO 6 J=1,NN 
DO 6 K=1,IBW 
BN(JU,K) = CN(J,K) 
SN(JU,K) = CN(JU,K) + DTIME*SN(J,K) 
WRITE (IOUT,607) 
RETURN 


Dh KK oR Ko ok OK KK ok OK KK ok ok ok 2 ok OK OK OK OK OK OK KOK KK KOK OK ok ok OK OK OK ok ok KK ok ok OK OK OO KOK KOK ROK KOK KO OK OK OK KO 


FORMAT STATEMENTS 

2K OK OK OK OK OK OK OK KK OK ok Kk ok ok ok KK KK OK OK OK OK OK KK OF OK OK OK OK KO OO OK OK kK OK OK OK OK OK OK KOK OK OF KOK OF OK OK OF OK OK OK KOK OK OF 

FORMAT (F10.0) 

FORMAT (15) 

FORMAT (1616) 

FORMAT (’0’,15X,’TIME INCREMENT =’,E13.6) 

FORMAT (’0’,15X,’NUMBER OF OUTPUT TIMES=’,15) 

FORMAT (’0’,15X, ‘MULTIPLES OF TIME INCREMENT AT WHICH RESULTS ‘’, 
‘ARE BEING PRINTED’ ) 

FORMAT (’ ’,15X,1616) 

FORMAT (’0’,15X,’EULERS METHOD IS BEING USED’) 

FORMAT (’0’,15X,’CRANK-NICHOLSON METHOD IS BEING USED’) 

FORMAT (’0O’,15X,’PURE-IMPLICIT METHOD IS BEING USED’) 

END 
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SUBROUTINE MODSN 
IMPLICIT REAL*8(A-H,0-Z) 

3K OK OK KK ok Ok OK ok ok 2k OK ok OK ok OK ok KOK ok ok ok ok ok ok ok ok ok Ok ok 
* READ SPECIFIED TEMPERATURES * 
* AND MODIFY MATRICES * 


2K OK KK KK Ko KKK KK ok OK ok ok OK OK KOK KO OK OK 


Dk OK KK RK OK KK kK ok KK ok ok ok ok kK ok Kok ok 2k ok ok OK ok ok Rk ok ook ok ok oR ok ok ok oko ok ok ok ok ok ok kok ok ok ok ok oF ok ok OK 


DIMENSIONS 

eK OK KOK ok OK ok ok ok Ok oko 2K ok 2K ok ok ok ok ok ok ok Ko OK ok ok kK ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok oF ok ok ok ok 
COMMON BN(500,30),CN(500.30),SN(500,30),SP(500),Y(500),xX(500), 
*T(500) ,GN(500),RN(500) , GEN(50),CEE(50),RHO(50) ,CAY(50),TIMAX, 
*DTIME, TIME .NTOUT(500) ,MN(500) ,NRL(500),KN(500),JUN(500),IN(500),uUC, 
“ ICOUNT ,NTIME,NT,NQ,NH,NC,NG,NE,NN,IBW,IOUT,IIN 


KOK KK OK EK KK KK OK Ko ok KK KR OK OK KK OF OK Ko ROK Ook ok oo OK KOR OK OK KOK KK OK Oo OF ok ob Ok Ok 


SPECIFIED NODAL TEMPERATURES 
KOK OK KOK OR OK OK OK OK OK OK OK OOK OK OF OK OK OOK OO OK OK OK OK OK OK KOK OOK OK OOK oO KO OK OOK OF OOO KOE KOK OO OK FOF Oh OK KOK OK OK 
WRITE (IOUT,605) 
DO 12 J=1,NN 
NRL(JU) = O 
DO 13 J=1,NT 
READ (IIN,504) K,T(K) 
NRL(K) = 1 
WRITE (IOUT,606) K,T(K) 


eK EO KOR ROR KOK RK KF oF RK oR RK eR Soi ok ak kok ok ok ok oF RR OR ok ko Kok ok bok ok kok 


MODIFY SN AND SP MATRICES FOR SPECIFIED NODAL TEMPERATURES 
DK OR OK eK oR OK Oo KK OK ok oF KK OK ok ok Ko Ok oR kK OK OK OK OK OK or OK Kk KF OK OO OK KOR ok ok ok ok oF Kok ok Ok ok Ok 
DO 14 J=1,NN 
IF (NRL(J)) 14,14,2 
We (USD) Wa %oe 
IL = 1 
LOS =a at 
th CU=TBW) 5e.5)4 
New RB Wee 
DOVGe = Te Ly 
Kae ite Abies st 
SP'(1) = SPCL)-=SNCL.K)*T Ca) 
SNCL KK) = 70,0 
ie CNN) aS, WEE 
KL = 2 
KU = IBW 
IF (J-NN+IBW-1) 10,10,9 
NU) SeNIN St alse 4] 
DO 11 K=KL,KU 


1 
I)-SN(J,K)*T(uJ) 
.O 

SINCUG a) & Ue 
CONTINUE 
RETURN 


SOK OR oR oR oR ak ok ok ok oko ao ok ok aK oR oo ok ok ok og ak ok ak aR Rok ok ok ok Rr kk ok OK oF ok ok kok 


FORMAT STATEMENTS 


SEK KORO ORR OKO OK OK IE ORR ROR I oR rR RR oR RR ok tok ok kk ok 
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504 FORMAT (15,F10.0) 

605 FORMAT (’0’,15X,’SPECIFIED NODAL TEMPERATURES’ / 
{ aS XC NOLES TEMP.’ ) 

606 FORMAT (7"", 15x13 F 10.3) 
END 
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SUBROUTINE DBAND(N,IBW,U,LABEL) 
IMPLICIT REAL*8(A-H,0-Z) 


DEE OK eK KK KK oo oo Ko Kok ok ok ok KOK KOK ook ok ok ok ok ok ok 


* DECOMPOSE SYMMETRIC BANDED MATRIX * 


EK EK OK EK KEE KK OK KK oo OK OK OK OK OK KK OK OK KK OK OK OK OK KOK kk 


3K OK OK KK KK OE KK KK ok ok OK OK KK OK KK OK OK OK OK ok oko KK ok KK OK KK KK ok ok ok ook ok ok ok ook oR oR oR KOK Ok OK ok ok Ook 


ARGUMENTS- 


N - NUMBER OF NODAL POINTS 
IBW - UPPER BAND WIDTH 
U - RESULT OF DECOMPOSITION ALGORITHM 


LABEL END OF FILE INDICATOR 
=O eR EXECURLON ES siORcONTENUE 


= th EXECUTION Sis) Os ERMINATE 


DEK OK KK A ok OE KK Ko Ko Kok KK ok KK EK KK Oo ok ko ok kK OK Kok OK a oo ok ok RO KF ROK OK ok ok oe ok ok 


DEK OK KE a ok KK OK OK OK kK ok KK OK KK XK ok ok ok ok ok ok ok OK OK Ko ok ok OK OK OK ok ob OK OK OK OF Ok KO ok OO OOK OK OK OK OK OK 


DIMENSION 


dK OK ok ok ok ok KK ok ok KK ok OK ok ok ok OK ok oh ok ok KK ok OK ok ko ok ok ok ok ok ok ok ok ok ok ok ok ok ok ok OK OK oF Ok OK KOK Ok OK OF OF OK Kk EOE 


DIMENSION U(500, 30) 


SKK OK ok ok ok ok oR ok oR oR ok ok ok KK ok ok ok ob OK ook ko ok ok kk ok kok kk KOK ROR RR EK RE EERE EKER EK 


DOVUEEE PREGCISTONSPSUMT Dis D2] D3s> DSORT 


DEK ok EO RK OK kK ok ok ok ok OOK Ko oo Ok ok oO KOK oko ok oo oR OK OR OK OK OR KK Kk ok ok Kk ok ok ok 


OK KK OK OK KK OK Oo KOK KOK EK OK oo ok ok ok ok OK ok RK ok OK ok ook ok OK KK oR OK OK Of ok OK RO Ok OF KO OK 


DECOMPOSE BAND (MAY REQUIRE DOUBLE PRECISION ACCUMULATORS. CHANGE 
SQRT TO DSQRT ON CARD FHT20649 AND ACTIVATE DOUBLE PRECISION 
DECLARATION ON CARD FHT20619. ) 
KK OK OK OK ok ok KO OK OK ok cK Ook ok ok Ko oo OK OK Ko ok ok OK kK OK OK ok OK OKO rok OK OO OK KO OO OF OK Ok Or OK OK 
DO 13 I=1,N 
es IN) Sah SP) 
IF (IBW - IP) 2,3,3 
IP = IBW 
bo 1 J=1,1P 
IQ = IBW - J 
ay XG tS te) Sa ae Boia 
TO ee 
DSUM = U(I,u) 
IF (10 - 1) 8,6,.6 
DOM7EK=1 10 
LMKe= Tok 
KP1 =K + 1 
Were Sf 2 tk 
Di = U(IMK,KP1) 
D2 = U(IMK,JPK) 
DSUM = DSUM - D1*D2 
me Cl Sa Se aie 
IF (DSUM) 10,10, 11 
LABEL = 1 
GOmnG 4 
D3 = 1.0/DSQRT(DSUM) 
U(I,JU) = D3 
GO TO 13 
U(I,J) = DSUM*D3 
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CONTINUE 
EABEE “= 
RETURN 
END 
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SUBROUTINE RHS 
PIMBE LCI REAE*SiCA-H OZ.) 


3 kok KOK Koo oko oR Ko Ko de ok ot ok ok ok oF ok 


* UPDATE RIGHT HAND SIDE * 


OK ROK RK RR RK ok RO ok oF ok ok ok ok ok Ok 


FOR KK OK oo Ko Ko KEK ok KK RK EK OK Koo Ko KOK ok ok oR oo KK KR RK KK OK eK oR ok oF ok ok ok ok 


DIMENSIONS 

2K OK KK OK OK OK ok 2K ok ok kK Kok OK OK KK OK OK ok KK oo 2K OK ok ok ok ok KK OK OK OK ok ok ck OK ok cK of ok oF OK ok KF OK ok OK OF OK ok OF OK ok ok KO Ok 
COMMON BN(500,30),CN(500, 30),SN(500,30),SP(500),Y(500),xX(500), 
*T(500),GN(500) ,RN(500) ,GEN(50),CEE(50),RHO(50) ,CAY(50),TIMAX, 
*DTIME,TIME,NTOUT(500) ,MN(500) ,NRL(500) ,KN(500),UN(500),IN(500),uUC, 
* ICOUNT ,NTIME,NT,NQ,NH,NC,NG.NE,NN,IBW,IOUT,IIN 


SKK KE KK OK KK KK ok KK oO OK ok kok ok KK oR Ok oo ok ok bok oR ok cK Oo KK KK OK OK OK OK OK OK KR OK ok OK bok 


STEADY STATE PROBLEMS 
POOR OR ROR RG UE OK ok GR I I aE Gk aE Rt ak ok ke de ak 
il? UNG)” AO), Los Be 
DO 21 J=1,NN 
RN(J) = GN(J) + SP(J) 
GORTORS 


OKO KK KOK KO OK KK ok KK OK KK oF Ko KK KK ok ok OK OK Ko ok ok OK oR oR KK ok ok ok ok aK kK oR KE KK XE ok oR ok ok ok ok 


INCREASE TIME 

FEO ER OR oR RR oR RE RI SER kok RR oR EER RRR Kok oR RR tok kK ER oh RR kok ok RRR RK + 
TIME = TIME + DTIME 

ICOUNT = ICOUNT + 1 


dk kK ok ok ok ok ok ok kK ok oo oko KOR OK OK ok Ok OK OR KO OK kK or Kk Ok oh ok ok ok oo KE Ok KK KK EO oR KK OK OF OK OF OE 


COMPUTE NEW RIGHT HAND SIDE 
Se OK ER OR RR OR RE OR oR RRR oR RR ok oe 
DO 8 J=1,NN 
RN(J) = GN(JU)*DTIME + SP(UJ) 
MES 4 
KU = IBW 
IF (J-NN+IBW-1) 2,2,1 
KU =a INNO 
DO 3 K=KL,KU 
Ulu) es dh se Ss 4 


RN(J) = RN(J) + BN(JU,K)*T(JJ) 
IF (JU-1) 8,8,4 
ML a & 
KU = IBW 
IF (wv = IBW) 5.6.6 
KUR=ad 


DO 7 K=KL,KU 
Gule = ee Kes 
RN(J) = RN(J) + BN(JuU,K)*T(UJ) 
CONTINUE 
RETURN 
END 
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SUBROUTINE MODRN 

IMPLICIT REAL*8(A-H,0-Z) 

KKK KO KK OK KKK KKK KK KKK KOK KE KOK KKK 
* MODIFY RIGHT HAND SIDE FOR * 
* SPECIFIED TEMPERATURES sg 


2K OK KK OK OOK a KK Oo Oo Kk OKO OK OF OK OOK KOK OK OK OK 


OK KK KOK OK ok ok OK OK oR KK ok Ko eo oR oo OK KK OK OK OK KK OK OK OO KK OO KOK 


DIMENSIONS 

KO OR KK OK OK KK OK KK ok KR or eK oO Kok ok ok ok ook ok oO Ko oR ok ok ok ook OR OK KK oR ok ok 
COMMON BN(500,30),CN(500,30),SN(500,30),SP(500),Y(500),xX(500), 
*T(500),GN(500),RN(500),GEN(50) ,CEE(50) ,RHO(50),CAY(50),TIMAX, 
*DTIME, TIME ,NTOUT(500),MN(500) ,NRL(500) ,KN(500),UN(500),IN(500),UC, 
be ICOUNT ,NTIME,NT,NQ,NH,NC,NG,NE,NN,IBW,IOUT,IIN 


OK KK KK KK KK KK OK oO KK OK KO RK ok KK ok OK ok ok ok ok ok oR ok oR OK OK ok OK oR OK KK OK OK OF ok ok Ok ok ok Ok oF 


MODIFY RN MATRIX FOR SPECIFIED NODAL TEMPERATURES 
Se OK ROK OK OK KO ORK OK ROR KOK KOKO RK KK KK OK KK oo ok OK ok ok ok ok oh ok OR RK OK oh OK OF ok OF 
DO 2 J=1,NN 

IF (NRECI)) 2.2.4 

RN(JU) = T(J) 

CONTINUE 
RETURN 
END 
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SUBROUTINE SBAND(N,IBW,U,B,X) 
IMPLICIT REAL*8(A-H,0-Z) 


SKK OK Ko aK ook ok ok ok kk kok ok ok Ok KOK KK 


* SOLVE SYMMETRIC BANDED MATRIX * 


OK eK KE OK OK KKK KK ok ok Kok ok ok KK ok ok ok ok 


3K KOK ok ok KK aK KK KK ok Ko ok OK ok ok ok ok oko oo ok ok kok kok ok kok ok ok ok ok ok ok ob Fok oo ok kok 


ARGUMENTS- 


N - NUMBER OF NODAL POINTS 
IBW - UPPER BAND WIDTH 


U - VALUE OF CONDUCTION MATRIX ELEMENT 
B - VALUE OF RIGHT HAND SIDE ELEMENT 
X - RESULT OF SOLUTION ALGORITHM 


SRK KK EO RK KK KK oR KK KK ok ook ok ok a ok ok ok kk kk ok oo OK KR KR RK RK ok Ook OK Kok ok ok 


DEK KK kK OK a OK aK KKK RK RK ROK OK KR ok oR ok of ok Ro ok Kook RoR Ook ok ok ok ok ok ok oo ok ok ok oe ok ok 


DIMENSIONS 


OK EK KK RK OK KK KK KK oo Ko ok ok oe KK OK ok ok oR ok ok ok ok ok ok ok ok ook ok ok ok ok ook ok ok ok ok ok ok ok ok ok 


DIMENSION U(500,30), B(500), xX(500) 


OK KK KR ok ok ok Oo EK OK oo ko oe KKK KK EK KK KR ok ok kok ok oo ok of ok ok ook oF ok ok oR ok 


DOVERE PREGES TON DSUMEND i e2 


SEK EK EK RK oR KR OR RR OR KOK oF oF ok oR ok oR ok ok bok ok ok ob ok ok oF KK OK Rob KF ROR EO 


2K KOK OK OK OR KK KK Kk ok KO OR OF ok oo kK ok eK ok KK ok OK ok ok Ook OK OF OK ok ok ok oF ok ok ok OK OK OK OK OE OK KR RO KKK KK EES 


SOLVE BAND (MAY REQUIRE DOUBLE PRECISION ACCUMULATORS. ACTIVATE 
DOUBLE PRECISION DECLARATION ON CARD FHT20762. ) 
2K OK OK OO OK OK OK OK OK OK OK Ook Oo KK OK OK OK OK OKO OF OK OO OK OK KOK KOK Ook OKO OKO OO OK OOO OK OO OO 
DO .5,.T=4,N 
Ses T -— LBW 4+ 94 
Le (Lt i= ew) 46.4.2 
Vesey 
DSUM = B(I) 
Kt =e } = 4 
IF (JU = K1) 3.3.5 
DO 4 K=J,K1 
LMKP1.= 1 = Kt 4 
Di = U(K,IMKP1) 
DSUM = DSUM - D1i*X(K) 


X(I) = DSUM*U(I,1) 
DO 10 L1=4 ,N 

bateNclter a 

J= I+ IBW - 1 

Tee CU=N) Vaan 

J o=aN 
DSUM .= eX-Gh) 
KD ees 


KMIP1 =K- I+ 1 

D2 = U(I,KMIP1) 

DSUM = DSUM - D2*X(K) 
X(I) = DSUM*U(I,1) 
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SUBROUTINE RESUL(IFLAG) 
IMPLICIT REAL*8(A-H.0-Z) 


2K KK OK OK KK ok KK OK Ok OK OK 


OP RUN RESULTS s* 


de OK OK KK OK KOK ok Ok OKO ok Ok 


3K OK EK OK ok Ko KK KE KK KK ok KK ok ok ok ok OK oF ok OK KR oko Kok ok ok kok ko kk ok ok ok ok ok kok ok ok ok ok ok 


ARGUMENT - 


IFLAG - PRINTOUT TIME FLAG 

OLE OULEUT ST LMESSPECTELCATLONSSARESGORREG Ti 

1 IF CALCULATIONS EXTEND BEYOND THE RANGE OF THE OUTPUT 

2 IF THE NUMBER OF OUTPUT TIMES DGES NOT MATCH THE 
NUMBER OF MULTIPLES OF THE TIME INCREMENT INPUT 


OK KK KO RK KK KO KK OK OK OOK OK oh ok oo KOKO OK OK KO KE KKK KK KOK KK KY 


KO KK KK KK OK OK OK KE OK KK oo ok ok ok ok ok ook ok ok Ok ok ok Ko OR oR OK OR OR OK OK Ok OOF 


DIMENSIONS 

FR RR ROKK ROE ERR RO OR OR OR IO Rf RR OOF a dak kk de okook kk a ok aK ear ok ok 
COMMON BN(500,30),CN(500,30),SN(500,30),SP(500),Y(500),X(500), 
*T(500),GN(500),RN(500) , GEN(50),CEE(50),RHO(50),CAY(50),TIMAX, 
*DTIME, TIME,NTOUT(500) ,MN(500) ,NRL( 500) ,KN(500) ,JUN(500),IN(500),uC, 
# ICOUNT ,NTIME,NT,NQ,NH,NC,NG,NE,NN,IBW,IOUT,IIN 


Sh OK OK ok ok ok ok ok KK 2K ok ok Kk OK ok ok OK OK ok OK OF OK Ok OF KOK KOR OK EK OKO KK KEKE EK KK EHR REE EK KEES 


STEADY STATE OR TRANSIENT PROBLEM 


KOK OK KK OK KK KK OK KK KOK KK KK OK oR KK OK OK oo ok oR OK OK OK KEK RK OK ok oR ok ok ok ok ob oR ok oF ok ok 


Le ANC} 443 


SKE ORO OOOO OR oR kok a Io a oo I RK oR ook or oR fo do ok ok ok ok ok ok kk ok ok 


PRINT NODAL TEMPERATURES - STEADY STATE PROBLEMS 
eR OF ROK ORO OR OR ORO OOO ROO RO ORR RRO RO ROR EEF ok Rk dk oF bok ob kok eA 
WRITE (IOUT,601) 
DO 2 J=1,NN 
WRITE (IOUT,602) JU,T(u) 
GO TO 7 


Sh KK ook oo ko oR ook oF ok ook ok KK Ko FR ok ok ok oR ok oR OK OO KK KK OK OK OR OF OK Ok OF kOe 


PRINT NODAL TEMPERATURES - TRANSIENT PROBLEMS 
sk ok OK OK OK ok oF ok kok ok ok ok ok ok ok ok ok ok ok or ok ok ok ok ao a ok ok ok ok ok of ok ok oR oo ok oR ok ok ok oR oF Oo ok or ok 
IF (TIME - DTIME) 4,4,5 
WRITE (IOUT,.603) 
CONTINUE 
IF (ICOUNT - NTOUT(JC)) 7,10,9 
WRITE (IOUT,604) TIME 
DO 6 J=1,NN 
WRITE (IOUT,605) J,T(J) 
JG = JC°+ 1 
TF (JC =<-NTUME = 47. toe 14 
CONTINUE 
RETURN 


SK ok OK OK KK OK Ko eo RK oO OK Kok ok ok oo oR eo ok ok OK OF OR KO Ok OK 


SET PRINTOUT TIME FLAG 
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* * 
* Program ‘'‘Vartest5' * 
* * 


* * * *€ &© *€ *& * € * ER EK SK 


This program calculates the variables required 
as input for the finite element program which 
couples heat and fluid flow in a heated oil 
sand foundation. 


The variables calculated by this program are: 


(i) Total Head (H) 

(ii) Hydraulic conductivity (k) 
(iii) Specific storage coefficient (Ss) 
(iv) Generation rate (g') 
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* MAIN PROGRAM * 
Ciotit Stell Sulet Steet Sa ae a eee ee ee 


Dimensioning of arrays 


REAL T(500),AW(500) ,AB(500) ,AST(500) ,AS(500) 

REAL MW(500),MWSPC,MB(500) ,DT(500) ,DU(500) ,MV 

REAL GAMMAW(500),GAMMAB(500) ,GAMMAF (500) , TPREV(500) 

REAL SS(500),G(500) , VISCO2(500) , VISCH4(500) ,vSITU(500) 
DIMENSION IT(500),N(500) ,JPR(500) ,I1D(500) ,x(500) ,¥(500) 
DIMENSION 2Z(500),KN(500) ,THEAD(500) ,NEG(500) ,NGL(500) 
DIMENSION NTE(500) ,NCUR(500) ,FAC(500) ,STRT(500) ,KNN(500) 
DIMENSION PERM(500) ,HED(20) ,M(500),NN(500) ,WAT(500),BIT(500) 
DIMENSION VOW(500) , VOB(500) , VOWW(500) , VOBB(500) , vISW(500) 
REAL NW,NB,MMW(500) ,MMB(500),VISF(500) , PPERM(500) ,SSEL(500) 
DIMENSION NODE(500) ,DPHEAD(500) , PERMEL(500) ,GEL(500) 

REAL VOBVS(500),BITVS(500) , VOWVS(500) ,WATVS (500) ,GABVS(500) 
REAL VOBBVS(500) , VOWWVS(500) ,DGAMF (500) ,GAFVS(500) ,GAWVS(500) 
REAL NS,D(500) , TSTRES (500) ,DGAMB(500) ,DGAMW(500) , THIN(500) 
REAL TNUT(500) ,NODEN (500) 


Set Material Constants 


GAMMAS= 18.0 


. TIMEIN=0.264E+07 


Read elevation heads for each node from 
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file #2 - "“ELEV.HEAD2" and Calculate 
Depth to node 


DO 3000 I=1,240 
READ(2,3001)2Z(1) 
D(I)=100.0-2(1) 
FORMAT(G20.0) 


Read New Formation Temperatures 
from file #3 - "TEMP.IN2" 


a a a we eee eee ee 


DO 3002 I=1,260 
READ(3,3003)NODEN(I) ,TNUT(I) 
IF (TNUT(I) .LT.5.0)TNUT(I)=5.0 
FORMAT (2G20.0) 


LSOD=- 1 
LOG=0 
DO 3411 J=1,20 
LSOD=LSOD+ 1 
DO 3412 Le 1 712 
LOG=LOG+ 1 
LSOD=LSOD+ 1 
NODE (LOG ) =NODEN (LSOD) 
T (LOG) =TNUT (LSOD) 
CONTINUE 
CONTINUE 


Read Previous Formation Temperatures 
from file #4 - "TEMP.PREV2" 


DO 3004 1=1,240 
READ(4,3003)NODE(I1),TPREV(I) 
IF(TPREV(I).LT.5.0)TPREV(I)=5.0 
REWIND 4 


Read New Total Heads From 
File #5 - "TOTAL.HEAD2" 


DO 6000 I=1,240 
READ(5,6001)THEAD(I ) 
FORMAT (G20.0) 


Store New Formation Temperatures 
in file #4 - "“TEMP.PREV2" 


DO 3005 I=1,240 
WRITE(4,3006)NODE(I),T(I) 
FORMAT(I5, |. ;f 10205200 9 


Variable Calculations 


DO 1 I=1,240 
TEMP=T (I) 
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TEMPVS=TPREV (I ) 


Calculate change in temperature at each node 


DT(1I)=T(1I)-TPREV(I) 


Calculate Incremental coefficient of 
Thermal Expansion for Water 


CALL ALPHAW(TEMP,AWSPC) 
AW(I)=AWSPC 


Calculate Incremental Coeffiecient of 
Thermal Expansion for Bitumen 


ew a ae ae a i a a a a i a ew ee ee ee ee 


CALL ALPHAB(TEMP,ABSPC) 
AB(I)=ABSPC 


Calculate Incremental Coeffiecient of 
Thermal Expansion for Soil Structure 


CALL ALPHST(TEMP,ASTSPC) 
AST(I )=ASTSPC 


Calculate Incremental Coefficient of 
Thermal Expansion for Sand Grains 


CALL ALPHAS(TEMP,ASSPC) 
AS(I)=ASSPC 


Calculate Coefficient of Compressibility 


of Water 


CALL MWATER(TEMP,MWSPC) 
MW (I) =MWSPC 


Calculate Cummulative Volume Change of 
Bitumen at new Temperature 


CALL TEXBIT(TEMP,BITSPC) 
BIT(1)=BITSPC 


Calculate Cummulative Volume Change of 
Bitumen at Previous Temperature 


CALL TEXBIT(TEMPVS,BITS) 
BITVS(1I)=BITS 


Calculate Cummulative Volume Change of 
Water at new Temperature 
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CALL TEXWAT (TEMP, WATSPC) 
WAT (I )=WATSPC 


Calculate Cummulative Volume Change of 
Water at Previous Temperature 


ee ww we we we ee we ew ew ew ew wm wm wm ew ew em wm em wm wm Om ew eH = 


CALL TEXWAT(TEMPVS, WATS ) 
WATVS (1) =WATS 


Calculate Unit Weight of Bitumen and 
Water at New Temperature 


a iw a i we a i a ww ww we 


VOB(I)=100.+BIT(I) 

VOW(I)=100.+WAT(I ) 

GAMMAW (I )=981./VOW(I ) 

GAMMAB(I)=1010.43/VOB(I ) 

VOBB(1I)=VOB(I)#*0.886 

VOWW(I)=VOW(1I)#0.114 

GAMMAF (I )=VOBB(I) /(VOBB(I )+VOWW(I) ) *GAMMAB(I)+ 
VOWW (I) /(VOBB(I )+VOWW(I) ) *+GAMMAW(I ) 


Calculate Unit Weight of Bitumen and 
Water at Previous Temperature 


ee ee we wwe a a a iw a i a a a ai a ee ee eS 


VOBVS (I )=100.+BITVS(I) 

VOWVS (1 )=100.+WATVS (I) 

GAWVS (I )=981./VOWVS (I) 

GABVS (I )=1010.43/VOBVS (I) 

VOBBVS(I1)=VOBVS(I1)*0.886 

VOWWVS (I )=VOWVS(1I)#*0.114 

GAFVS (I )=VOBBVS (I) /(VOBBVS(I )+VOWWVS (I) )*GABVS(I)+ 
VOWWVS (1) /(VOBBVS (I) +VOWWVS (I ) ) *GAWVS (I) 


Calculate Change in Unit Weight due to 
Temperature Change 


DGAMF (I )=GAMMAF (I )-GAFVS (I) 
DGAMB (I )=GAMMAB(I )-GABVS (I) 
DGAMW (I )=GAMMAW(I )-GAWVS (I ) 


Calculate in situ Viscosity for Bitumen 


VISCH4(I)=(1.9661227*10.#*9, )*(T(I)**(-4.936318) ) 
VISCO2(1I)=(3.5066118*10.*##7.)*(T(1I)**(-4. 107852) ) 
VSITU(I )=(VISCH4(1I)#*.8)+(VISCO2(1I)*.2) 


Calculate Viscosity for Water 


CALL VISWAT (TEMP, VWSPC) 
VISW(I)=VWSPC 
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Calculate Viscosity for Fluid (bitumen & water) 


VISF(1I)=.886#VSITU(I)+.114#VISW(I) 


Calculate in situ Hydraulic Conductivity 


PERM (I) =(GAMMAF(1)#*3.0E-12)/(VISF(I)/1000.) 


Calculate Coefficient of Compressibility 
For Bitumen 


MB(I)=(-6.25E-9#T(I))+2.15E-6 


Calculate change in pore fluid pressure 


DU(I)=(NB*DT(I)#(AB(I)-AS(1I))+NW 
*DT(I)*(AW(1)-AS(1))-AST(1)#DT(I )+MV*(NB 
*DGAMB(I )#D(I )+NW*#DGAMW(I)#*D(1I)))/(MV+NB 
*(MB(I))+NW*(MW(I))) 


Set pore pressure change at surface nodes 
to zero 


IF(I.EQ.12)DU(I 
IF(I.EQ.24)DU(I 
IF(I.EQ.36)DU(I 
IF(I.EQ.48)DU(I 
IF(I.EQ.60)DU(I 
IF(I.EQ.72)DU(I 
IF(I.EQ.84)DU(I1 
IF(I.EQ.96)DU(I 
IF(I.EQ. 108)DU( 
IF CT SEO. 120). Dut 
( 
( 
( 
( 
( 
( 
( 
( 
( 
( 
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IF(I.EQ.132)DU 
IF(I.EQ.144)DU 
IF(I.EQ.156)DU 
IF(1.EQ.168)DU 
IF(I.EQ.180)DU 
IF(I.EQ.192)DU 
IF(1.EQ.204)DU 
IF(1.EQ.216)DU 
IF(I.EQ.228)DU 
IF(I.EQ.240)DU 


So. 2 Se SO Oo 
e e e e e 
DODDODOOCOOO0CO 


ee 


Calculate change in pressure Head 


ee 


DPHEAD(I )=DU(I ) /GAMMAF (I ) 


Calculate new total Head 


THEAD (I )=THEAD(I )+DPHEAD(I ) 
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Check if Pore Fluid Pressure > Total Stress 
if so, set Pore Fluid Pressure = Total Stress 
Calculate New Total Head 


we ww a ee a a we a a i a we we ee ee 


Total Stress due to tank 
14.63 m * 8.713 KN/m3 + .3048 m * 76.979 KN/m3 
150.93 KPa 


iow 


Total Stress due to Gravel Pad 
= 2m * 21.21 KN/m3 = 42.42 KPa 


ANANAANANAANANANANA 


TANK=0.0 

PAD=0.0 

IF(I.LE.96)TANK=150.93 

IF(I.LE. 156) PAD=42.42 

NS=1.0-NB-NW 

TSTRES (I )=D(I ) *(NB*GAMMAB (I ) +NW*GAMMAW(I )+NS*GAMMAS )+TANK 
& +PAD 

IF (((THEAD(I)-Z(1))*#GAMMAF(I)).LE.TSTRES(I))GO TO 9000 

THEAD(I)=TSTRES (I) /GAMMAF(I)+Z(I) 


Calculate generation rate 


ee a a i ww a a a ew ee ee 


ANAND 


9000 G(1)=(NW*AW(I) #DT(1)+NB*AB(I)#DT(I) 
& +AS(1I)#*DT(1I)#*(NB+NW)-AST(1I)*DT(1I))/TIMEIN 


Calculate specific storage coefficient 


NAAN 


SS(1)=( (NW*MW(I ) *GAMMAW(I ) )+(NB*MB(I ) *GAMMAB(TI ) ) 
& - (MV*GAMMAF (I ) ) ) *GAMMAF (I ) 
1 CONTINUE 


Set Generation rate for surface nodes equal to 
zero 


AQaGG a 


DO 51 J=12,240,12 
51 G(J)=0.0 


specific storage for each element 


c 
Cc Calculate hydraulic conductivity, generation rate and 
Cc 
(e 


WRITE(9,3085) 
3085 FORMAT(3X,'ELEM',2X,'NODE',2X,'NODE',2X,'NODE', 2X, 
*'NODE',2X,'PERMEL', 10X,'GEL',12X,"SSEL',//) 


DO 4000 I=1,209 
READ(7,4001)MNK1,MNK2,MNK3,MNK4 
PERMEL (I ) = (PERM(MNK1)+PERM(MNK2)+PERM(MNK3) 
& +PERM(MNK4))/4. 
GEL (I )=(G(MNK1)+G(MNK2)+G(MNK3)+G(MNK4) )/4. 
SSEL(I )=(SS(MNK1)+SS(MNK2)+SS(MNK3)+SS(MNK4))/4. 
WRITE(9,2298)1,MNK1,MNK2,MNK3,MNK4, PERMEL(I) ,GEL(I) ,SSEL(I) 
2298 FORMAT(3X,515,3X%,E12.6,3X,E12.6,3X,E12.6) 
4000 CONTINUE 
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FORMAT (415) 


PRINT OUTPUT 


DATA CHECK 


WRI TE(9,9600) 

FORMAT (Bitar ox T" Bx ST! 0X TPREV 9X. DT’. 13k, .D” 14x, 
“UAB ¥AZee AS' 7 12k, 'AW', 10%, "AST ,//) 

DO 8989 I=1,240 
WRITE(9,96)1,T(1),TPREV(I),DT(I),D(I),AB(1),AS(I),AW(I), 
&AST(I) 

CONTINUE 

WRITE(9,9700) 

FORMAT('1',2X,'I',5X,'DGAMB' ,9X,'DGAMW',12X,'MB',10X,'MW', 
*12X,'DU',12X, 'GAMMAF' ,//) 

DO 8976 I=1,240 

WRITE(9,96)1I ,DGAMB(I) ,DGAMW(I),MB(I),MW(I),DU(I) ,GAMMAF(I) 
WRITE(9,9800) 
FORMAT('1',2X,'I',6X,'DPHEAD',8X,'TSTRES' ,8X,'THEAD', 10X, 
4G UFI2KESSS® 12k.’ PERM" //) 

DO 8977 I=1,240 
WRITE(9,96)1,DPHEAD(I),TSTRES(I),THEAD(I),G(1I),SS(1),PERM(I) 
FORMAT(GS5,2X,E12.6,2X,E12.6,2X,512.6,2X,E12.6,2X,E12.6, 
F2KFE 12.6: 2%,612 6,2%,E12<6) 


READ(1,100) (HED(I) ,IJ=1, 16) 

FORMAT (20A4 ) 
WRITE(6,100)(HED(1),I=1,16) 
READ(1,101)NUMNP,NEGL,NEGNL,MODEX,NPER,TSTART,IPRI ,ITP56 
FORMAT(515,F10.2,215) 
WRITE(6,101)NUMNP,NEGL,NEGNL,MODEX,NPER,TSTART,IPRI,ITP56 
READ(1,102)IHEAT, IHEATN 

FORMAT (215) 

WRITE(6,102)IHEAT, IHEATN 
READ(1,103)IEIG 

FORMAT (15) 

WRITE(6,103)IEIG 
READ(1,104)ISREF,IEQUIT,ITEMAX,RTOL 
FORMAT (315,E10.4) 
WRITE(6,104)1ISREF,IEQUIT,ITEMAX,RTOL 
READ(1,105)IOPE,ALPHA 
FORMAT(110,F10.1) 
WRITE(6,105)IOPE,ALPHA 
READ(1,106)NPB 

FORMAT(15) 

WRITE(6,106)NPB 

READ( 1,107) IPNODE 

FORMAT(15) 

WRITE(6, 107) IPNODE 

READ(1,108)NSPER 

FORMAT (15) 

WRITE(6,108)NSPER 

READ(1,109)DTPER 

FORMAT(F 10.2) 

WRITE(6,109)DTPER 

po 1000 I=1,240 
READC1,770)19T(1) NGI), JPRCI ),IDCI):, 201), 2(1) -2(1) ,KN(T) 
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FORMAT(A1,14,A1,14,3F10.1,15) 
WRITE(6,110)IT(I),N(1),JPR(I),ID(I),X(1),Y(1),Z(1),KN(I) 
CONTINUE 
READ(1,111)ICON,IPRIC 
FORMAT (215) 
WRITE(6,111)ICON,IPRIC 
LNNN=0 
LMNO=- 1 
DO 9990 J=1,19 

LMNO=LMNO+ 1 

DO 9991 I=1,11 

LNNN=LNNN+ 1 

LMNO=LMNO+ 1 

THIN (LNNN ) =THEAD(LMNO) 
CONTINUE 
CONTINUE 
WRITE(6,112) (THEAD(K) ,K=1,240) 
FORMAT (6E12.6) 
READ(1,113)NTTF,NPTM, ILCOV,NTEMP,NOCV,NORA,NLOAD,NIHT,NTOT 
FORMAT (915) 
WRITE(6,113)NTTF,NPTM, ILCOV,NTEMP,NOCV,NORA,NLOAD,NIHT,NTOT 
DO 9944 KK=1,2 
READ(1,114)NTF,NPTS 
FORMAT (215) 
WRITE(6,114)NTF,NPTS 
READ(1,115)TIMV1,RV1,TIMV2,RV2 
FORMAT (4F10.1) 
WRITE(6,115)TIMV1,RV1,TIMV2,RV2 
DO 4465 IB=1,2 
READ(1,4466)NUD,NCOR,FUC, BSTRT,KNKN 
FORMAT(215,2F10.1,15) 
WRITE(6,4466)NUD,NCOR,FUC, BSTRT,KNKN 
CONTINUE 
READ(1,116)NEG(I),NGL(I),NTE(I) 
WRITE(6,116)NEG(1),NGL(I),NTE(I) 

DO 1003 J=1,209 
READ(1,117)M(J),NCUR(J),FAC(J),STRT(J),KNN(J) 
IF(GEL(J).EQ.0.0)NCUR(J)=2 
WRITE(6,127)M(J),NCUR(J),GEL(J),STRT(J) ,KNN(J) 

CONTINUE 

FORMAT (315) 

FORMAT (215,2F10.0,15) 

FORMAT(215,E10.4,F10.1,15) 
READ(1,118)NPAR1,NPAR2,NPAR3,NPAR4,NPAR5,NPAR7,NPARIO, 


&NPAR15,NPAR16,NPAR17,NPAR18 


FORMAT (514,18,112,120,314) 
LLL=0 

KKK=0 

DO 2000 J=1,2 

NNNN=209 
IF(J.EQ.2)NPAR2=19 
IF(J.EQ.2)GO TO 2000 
IF(J.EQ.2)NPAR16=1 
IF(J.EQ.2)NPAR3=1 
IF(J.EQ.2)NPAR4=1 
WRITE(6,118)NPAR1,NPAR2,NPAR3,NPAR4,NPAR5,NPAR7,NPARIO, 


&NPAK15,NPAR16,NPAR17,NPAR18 


DO 1008 I=1,NNNN 
LLL=LLL+1 
IF(LLL.GT.210)GO TO 1008 
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READ(1,119)NN(LLL) 

CONTINUE 

DO 1007 L=1,NNNN 

KKK=KKK+1 

IF(KKK.GT.210)GO TO 1007 

WRITE(6,119)NN(KKK) 

FORMAT (15) 

WRITE(6, 122) PERMEL (KKK) 

FORMAT(E10.4) 

WRITE(6,122)SSEL(KKK) 

CONTINUE 

DO 1005 K=1,NNNN 
READ(1,120)M2,IEL,MTYP,KG,BET,THIC,ETIME 

FORMAT (415, 3G10.0) 
WRITE(6,128)M2,IEL,MTYP,KG, BET, THIC,ETIME 

FORMAT (415,3F10.1) 
READ(1,121)NOD1,NOD2,NOD3,NOD4,NOD5,NOD6,NOD7,NOD8 
FORMAT (815) 
WRITE(6,121)NOD1,NOD2,NOD3,NOD4,NOD5,NOD6,NOD7 , NOD8 
CONTINUE 

CONTINUE 

WRITE(6,10) 
FORMAT(12X,'T(C)',6X, 'ALPHAW' ,9X,'ALPHAB',9X'ALPHAST', 
*8X'ALPHAS' ,9X'MWATER' ) 

DO 2 1=6,200 
WRITE(6,11)T(1),AW(I),AB(I),AST(I),AS(1I),MW(I) 
FORMAT(10X,F5.1,4F15.8,F16.9) 
READ(1,7777)NUT1,NUT2,NUT3,NUT4 ,NUT5,NUT6,NUT7,NUTS8,NUTS, 
&NUT 10 

FORMAT(14,18,414,18,20X,14,214) 
WRITE(6,7777)NUT1,NUT2,NUT3,NUT4,NUT5,NUT6,NUT7,NUT8,NUTO, 
&NUT 10 

READ(1,7778)NICE,CONVN 

FORMAT(15,E10.4) 

WRITE(6,7778)NICE,CONVN 

DOM T 9 ri t=4 2 
READ(1,7780)MUD,NODC1,NODC2,NODC3,ICTYP,KGLL,THUC, EATIME 
FORMAT(615,2F10.1) 
WRITE(6,7780)MUD,NODC1,NODC2,NODC3,ICTYP,KGLL,THUC,EATIME 
CONTINUE 

PRINT 2001 

FORMAT ('STOP' ) 

STOP 

END 


* * * * * * * *€ KO KOE EK OK OF 


* SUBROUTINE ALPHAW * 
* * * * & KK K€ KK K EK K 


This subroutine calculates the incremental 
coefficient of thermal expansion of water 
for a given temperature 


SUBROUTINE ALPHAW (TSPEC,AWSPC) 

DIMENSION AW(35),TW(35) 

DATA TW(1)/5.5/,TW(2)/8./,TW(3)/9.5/,TW(4)/10.5/, 
#TW(5)/11.5/,TW(6)/12.5/,TW(7)/13.5/,TW(8)/14.5/, 
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#TW(9)/15.5/,TW(10)/16.5/, mit 11)/17.5/,TW(12)/18.5/, 
#TW(13)/19. 5/5 TW(14)/21./,TW(15)/23.5/,TW(16)/26./, 
*#TW(17)/28.5/,TW(18)/35./,TW(19)/45./,TW(20)/55./, 
*TW(21)/65./,TW(22)/75./,TW(23)/85./,TW(24)/95./, 
*TW(25)/105./,TW(26)/115./,TW(27)/125./, TW(28)/135./, 
*TW(29)/145./,TW(30)/155./,TW(31)/165./,TW(32)/175./, 
*TW(33)/185./,TW(34)/195./,TW(35)/205.1/,AW(1)/3.34E-5/, 
#AW(2)/5.E-5/,AW(3)/1.001E-4/,AW(4)/1.001B-4/,AW(5)/1.E-4/, 
*AW(6)/1.001E-4/,AW(7)/1.001E-4/,AW(8)/2.E-4/,AW(9)/2.002E-4/, 
*AW(10)/2.003E-4/,AW(11)/1.993E-4/,AW(12)/2.002E-4/, 
*AW(12)/2.002E-4/,AW(13)/2.001E-4/,AW(14)/2.496E-4/, 
#AW(15)/2.663E-4/,AW(16)/2.493E-4/,AW(17)/2.66E-4/, 
*AW(18)/3.488E-4/,AW(19)/4.27E-4/,AW(20)/4.944E-4/, 
#AW(21)/5.609E-4/,AW(22)/6.164E-4/,AW(23)/6.71E-4/, 
*#AW(24)/7.341E-4/,AW(25)/7.768E-4/,AW(26) /8.278E-4/, 
*AW(27)/8.872E-4/,AW(28)/9.355E-4/,AW(29)/1.0009E-3/, 
*AW(30)/1.0553E-3/,AW(31)/1.1077E-3/,AW(32)/1.1855E-3/, 
*AW(33)/1.2603E-3/,AW(34)/1.3323E-3/,AW(35)/1.4274E-3/ 
AWSPC=-1.0 

IF(TSPEC-TW(1)) 90,90,91 

AWSPC=AW(1) 

RETURN 

IF(TSPEC-TW(35)) 93,92,92 

AWSPC=AW(35) 

RETURN 

DO 100"9=1 ,35 

IF(TSPEC-TW(J)) 100,94,100 

AWSPC=AW(J) 

RETURN 

CONTINUE 

MIN=1 

MAX=35 

CONTINUE 

IF((MAX-MIN)-1) 200,200,104 

NEXT= (MIN+MAX) /2 

IF(TSPEC-TW(NEXT)) 105,105,107 

MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN 

TWFAC=(TSPEC-TW(1I))/(TW(I+1)-TW(I) ) 
AWSPC=AW(I1)+TWFAC# (AW(I+1)-AW(I)) 

RETURN 

END 


* * * * * * * K€ KK KOK KOK OK OX 


* SUBROUTINE ALPHAB * 
* * * * € * * € # * EE KK 


This subroutine calculates the incremental 
coefficient of thermal expansion of bitumen 
for a given temperature 


SUBROUTINE ALPHAB (TSPEC,ABSPC) 

DIMENSION AB(19),TB(19) 

DATA TB(1)/15./,TB(2)/25./,TB(3)/35./,TB(4)/45./, 
#TB(5)/55./,TB(6)/65./,TB(7)/75./, 
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*TB(8)/85./,TB(9)/95./,TB(10)/105./,TB(11)/115./, 
SPBC120/125./, 7B \ 13)% 135 .7,7B (14) /145.7,TB015)/ 155.7, 
*TB(16)/165./,TB(17)/175./,TB(18)/185./,TB(19)/195./, 
*#AB(1)/1.0225E-3/,AB(2)/1.0185E-3/,AB(3)/1.0159E-3/, 
*AB(4)/1.0143E-3/,AB(5)/1.0142E-3/,AB(6)/1.0153E-3/, 
*AB(7)/1.0177E-3/,AB(8)/1.0212E-3/,AB(9)/1.0263E-3/, 
*#AB(10)/1.0323E-3/,AB(11)/1.0399E-3/,AB(12)/1.0486E-3/, 
*AB(13)/1.0585E-3/,AB(14)/1.0699E-3/,AB(15)/1.0824E-3/,. 
*AB(16)/1.0962E-3/,AB(17)/1.1112E-3/,AB(18)/1.1277E-3/, 
*AB(19)/1.1453E-3/ 

ABSPC=-1.0 

IF(TSPEC-TB(1)) 90,90,91 

ABSPC=AB(1) 

RETURN 

IF(TSPEC-TB(19)) 93,92,92 

ABSPC=AB( 19) 

RETURN 

DO 100 J=1,19 

IF(TSPEC-TB(J)) 100,94,100 

ABSPC=AB(J) 

RETURN 

CONTINUE 

MIN=1 

MAX=19 

CONTINUE 

IF((MAX-MIN)-1) 200,200,104 

NEXT=(MIN+MAX) /2 

IF(TSPEC-TB(NEXT)) 105,105,107 

MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN 

TBFAC=(TSPEC-TB(1I))/(TB(I+1)-TB(1) ) 
ABSPC=AB(I1)+TBFAC#(AB(I+1)-AB(I) ) 

RETURN 

END 
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* SUBROUTINE ALPHAST * 
* * * * * *€ € *€ € € € HK KE K 


This subroutine calculates the incremental 
coefficient of thermal expansion of soil 
structure for a given temperature 


SUBROUTINE ALPHST (TSPEC,ASTSPC) 

DIMENSION AST(14) ,TST(14) 

DATA -TST(1)/29.75/; TST(2)/23.757,TST(3)/56.25/, 
*TST(4)/68.75/,TST(5)/81.25/,TST(6)/93.75/, 
*TST(7)/106.25/,TST(8)/118.75/,TST(9)/131.25/, 
*TST(10)/143.75/,TST(11)/156.25/,TST(12)/168.75/, 
*TST(13)/181.25/,TST(14)/193.75/, 
*AST(1)/3.2E-5/,AST(2)/4.4E-5/,AST(3)/4.8E-5/, 
*AST(4)/4.4E-5/,AST(5)/4.4E-5/,AST(6)/5.2E-5/, 
*AST(7)/4.8E-5/,AST(8) /4.4E-5/,AST(9)/4.4E-5/, 
#AST(10)/4.0E-5/,AST(11)/4.0E-5/,AST(12)/4.8E-5/, 
*AST(13)/4.0E-5/,AST(14)/4.8E-5/ 
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ASTSPC=-1.0 

IF(TSPEC-TST(1)) 90,90,91 
ASTSPC=AST(1) 

RETURN 

IF(TSPEC-TST(14)) 93,92,92 
ASTSPC=AST(14) 

RETURN 

DO 100 J=1,14 

IF(TSPEC-TST(J)) 100,94,100 
ASTSPC=AST(J) 

RETURN 

CONTINUE 

MIN=1 

MAX=14 

CONTINUE 

IF ((MAX-MIN)-1) 200,200,104 
NEXT= (MIN+MAX) /2 

IF (TSPEC-TST(NEXT)) 105,105,107 
MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN 
TSTFAC=(TSPEC-TST(1I))/(TST(1+1)-TST(I) ) 
ASTSPC=AST(I)+TSTFAC# (AST(I1+1)-AST(1I) ) 
RETURN 

END 


* * * * *€ * * € * KE KK K 
% SUBROUTINE ALPHAS * 
* * * * * *€ * & * KK KE EK F 


This subroutine calculates the incremental 
coefficient of thermal expansion of sand 
grains for a given temperature 


SUBROUTINE ALPHAS (TSPEC,ASSPC) 

DIMENSION AS(5),TS(5) 

DATA TS6:1)/20.7 -TS(2)7 (oes te Gs) 7125.67 ,TS(4)/ 175.7, 
*TS(5)/225./, 
#AS(1)/3.4E-5/,AS(2)/3.8E-5/,AS(3)/4.E-5/, 
*#AS(4)/4.4E-5/,AS(5)/5.0E-5/ 

ASSPC=-1.0 

IF(TSPEC-TS(1)) 90,90,91 

ASSPC=AS (1) 

RETURN 

IF(TSPEC-TS(5)) 93,92,92 

ASSPC=AS(5) 

RETURN 

BO 100 J=1,5 

IF(TSPEC-TS(J)) 100,94,100 

ASSPC=AS (J) 

RETURN 

CONTINUE 

MIN=1 

MAX=5 

CONTINUE 

IF ((MAX-MIN)-1) 200,200,104 
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NEXT= (MIN+MAX) /2 
IF(TSPEC-TS(NEXT)) 105,105,107 
MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN 
TSFAC=(TSPEC-TS(1))/(TS(I+1)-TS(I)) 
ASSPC=AS (I )+TSFAC*(AS(I+1)-AS(I) ) 
RETURN 

END 


* * *¢+ © *€ *€ € € & EE Ee KE OF 


* SUBROUTINE MWATER * 
* * * *€ € * * € & ee eK EF 


This subroutine calculates the incremental 
coefficient of thermal compressibility of 
sand grains for a given temperature 


SUBROUTINE MWATER (TSPEC,MWSPC) 
REAL MW(5),TW(5) ,MWSPC 
DATA TW(1)/20./,TW(2)/100./,TW(3)/125./, 


*TW(4)/150./,TW(5)/200./, 
*#MW(1)/4.5E-7/,MW(2)/4.73E-7/,MW(3)/5.15E-7/, 
*MW(4)/6.04E-7/,MW(5)/8.54E-7/ 


MWSPC=-1.0 

IF(TSPEC-TW(1)) 90,90,91 
MWSPC=MW( 1) 

RETURN 

IF(TSPEC-TW(5)) 93,92,92 
MWSPC=MW(5) 

RETURN 

DO 100 J=#1,5 

IF(TSPEC-TW(J)) 100,94,100 
MWSPC=MW(J) 

RETURN 

CONTINUE 

MIN=1 

MAX=5 

CONTINUE 

IF((MAX-MIN)-1) 200,200,104 
NEXT=(MIN+MAX) /2 
IF(TSPEC-TW(NEXT)) 105,105,107 
MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN 
TWFAC=(TSPEC-TW(1I))/(TW(1I+1)-TW(I) ) 
MWSPC=MW(I)+TWFAC* (MW(1+1)-MW(I) ) 
RETURN 

END 


* * * * * * *€ € € EK KOK OK K OF 


* SUBROUTINE TEXBIT * 
* * * * * * &€ € * Ee KE KE EK F 
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This subroutine calculates the thermal 
expansion (volume change) of in situ 
bitumen for a given temperature 


SUBROUTINE TEXBIT (TSPEC,BITSPC) 

REAL BIT(9),TB(9),BITSPC 

DATA TB(1)/5./,TB(2)/22./,TB(3)/50./,TB(4)/75./, 
*TB(5)/100./,TB(6)/125./,TB(7)/150./, 
*TB(8)/175./,TB(9)/200./,BIT(1)/-1.4215/, 
BI T(2)/0.0/ sBIT(3) /2.651/7, BIT(4)/5.373/, 
*BIT(5)/7.975/,BIT(6)/10.514/,BIT(7)/13.158/, 
*BIT(8)/15.971/,BIT(9)/18.737/ 
BITSPC=-2.0 

IF(TSPEC-TB(1)) 90,90,91 

BITSPC=BIT(1) 

RETURN 

IF(TSPEC-TB(9)) 93,92,92 

BITSPC=BIT(9) 

RETURN 

DO 100 J=1,9 

IF(TSPEC-TB(J)) 100,94,100 

BITSPC=BIT(J) 

RETURN 

CONTINUE 

MIN=1 

MAX=9 

CONTINUE 

IF ((MAX-MIN)-1) 200,200,104 

NEXT= (MIN+MAX) /2 

IF(TSPEC-TB(NEXT)) 105,105,107 

MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN 
TBFAC=(TSPEC-TB(1I))/(TB(I+1)-TB(I) ) 
BITSPC=BIT(1)+TBFAC#(BIT(I+1)-BIT(I) ) 
RETURN 

END 


* * * * * * € KF KO KE K KE KOK HK OF 


* SUBROUTINE TEXWAT * 
* * * * € * KE KK KE KK KE K F 


This subroutine calculates the thermal 
expansion (volume change) of water 
for a given temperature 


SUBROUTINE TEXWAT (TSPEC,WATSPC) 

DIMENSION WAT(19),TW(19) 

DATA TW(1)/23./,TW(2)/30./,TW(3)/40./,TW(4)/50./, 
*TW(5)/60./,TW(6)/70./,TW(7)/80./, 
*TW(8)/90./,TW(9)/100./,TW(10)/110./,TW(11)/120./, 
*TW(12)/130./,TW(13)/140./,TW(14)/150./,TW(15)/160./, 
*TW(16)/170./,TW(17)/180./,TW(18)/190./,TW(19)/200./, 
*WAT(1)/0.0/,WAT(2)/0.182/,WAT(3)/0.531/, 
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*WAT(4)/0.960/,WAT(5)/1.470/,WAT(6)/2.029/, 
*WAT(7)/2.668/,WAT(8)/3.357/,WAT(9)/4.106/, 
*WAT(10)/4.924/,WAT(11)/5.793/,WAT(12)/6.731/, 
*WAT(13)/7.730/,WAT(14)/8.798/,WAT(15)/9.946/, 
*WAT(16)/11.174/,WAT(17)/12.492/,WAT(18)/13.900/, 
*#WAT(19)/15.428/ 

WATSPC=-1.0 

IF(TSPEC-TW(1)) 90,90,91 

WATSPC=WAT (1) 

RETURN 

IF(TSPEC-TW(19)) 93,92,92 

WATSPC=WAT( 19) 

RETURN 

DO 100 J=1,19 

IF(TSPEC-TW(J)) 100,94,100 

WATSPC=WAT (J) 

RETURN 

CONTINUE 

MIN=1 

MAX=19 

CONTINUE 

IF ((MAX-MIN)-1) 200,200,104 

NEXT= (MIN+MAX ) /2 

IF (TSPEC-TW(NEXT)) 105,105,107 

MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN ) 
TWFAC=(TSPEC-TW(1))/(TW(I+1)-TW(I)) 
WATSPC=WAT(I)+TWFAC# (WAT(I+1)-WAT(I) ) 
RETURN 

END 


* * * * * * EK KK KK KK K KK K OF 


* SUBROUTINE VISWAT * 
* * * * * * * € F KK KK KEK KE FE 


This subroutine calculates the viscosity 
of water for a given temperature 


SUBROUTINE VISWAT (TSPEC,VWSPC) 

DIMENSION VW(35),TW(35) 

DATA TW(1)/5./,TW(2)/10./,TW(3)/15./,TW(4)/20.0/, 
*TW(5)/25.0/,TW(6)/30.0/,TW(7)/35.0/,TW(8)/40.0/, 
*TW(9)/45.0/,TW(10)/50.0/,TW(11)/55.0/,TW(12)/60.0/, 
*TW(13)/65.0/,TW(14)/70./,TW(15)/75.0/,TW(16)/80./, 
*TW(17)/85.0/,TW(18)/90./,TW(19)/95./,TW(20)/100./, 
*TW(21)/105./,TW(22)/110./,TW(23)/115./,TW(24)/120./, 
*TW(25)/125./,TW(26)/130./,TW(27)/135./,TW(28)/140./, 
*TW(29)/145./,TW(30)/150./,TW(31)/160./,TW(32)/170./, 
*TW(33)/180./,TW(34)/190./,TW(35)/200./,VW(1)/1.501E-3/, 


*VW(2)/1.3E-3/,VW(3)/1.136E-3/,VW(4)/1.002E-3/,VW(5)/8.9E-4/, 
*VW(6)/7.97E-4/,VW(7)/7.18E-4/,VW(8) /6.51E-4/,VW(9)/5.94E-4/, 


#VW(10)/5.44E-4/,VW(11)/5.01E-4/,VW(12)/4.63E-4/, 
*VW(13)/4.30B-4/,VW(14)/4.E-4/,VW(15)/3.74E-4/, 
*VW(16)/3.51B-4/,VW(17)/3.3E-4/,VW(18)/3.11E-4/, 
*VW(19)/2.94E-4/,VW(20)/2.79E-4/,VW(21)/2.65E-4/, 
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* #VW(22)/2.52E-4/,VW(23)/2.41E-4/,VW(24)/2.3E-4/, 


90 


oI 
92 


93 
94 
100 


102 
104 
105 
107 
200 


*#VW(25)/2.2E-4/,VW(26)/2.11E-4/,VW(27)/2.03E-4/, 
*#VW(28)/1.95E-4/,VW(29)/1.88E-4/,VW(30)/1.81E-4/, 
*VW(31)/1.69E-4/, VW(32)/1.59E-4/,VW(33)/1.49E-4/, 
*#VW(34)/1.41E-4/,VW(35)/1.34E-4/ 
VWSPC=-1.0 

IF(TSPEC-TW(1)) 90,90,91 
VWSPC=VW(1) 

RETURN 

IF(TSPEC-TW(35)) 93,92,92 
VWSPC=VW(35) 

RETURN 

DO 100 J=1,35 

IF(TSPEC-TW(J)) 100,94,100 
VWSPC=VW(J) 

RETURN 

CONTINUE 

MIN=1 

MAX=35 

CONTINUE 

IF((MAX-MIN)-1) 200,200,104 

NEXT= (MIN+MAX) /2 

IF(TSPEC-TW(NEXT)) 105,105,107 
MAX=NEXT 

GO TO 102 

MIN=NEXT 

GO TO 102 

I=MIN ; 
TWFAC=(TSPEC-TW(1))/(TW(I+1)-TW(I) ) 
VWSPC=VW(1)+TWFAC# (VW(I+1)-VW(I) ) 
RETURN 

END 
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